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profile of a ship to the spatial frequency components. The 
shape of the coefficient curve can be used to classify the 
type of a ship from 8 different categories. Bigeye the 
differences is so minor that it is difficult to implememem 
computer proeqmkham  Gemueceoqnize 1:2. The second method is a 
B-spline Coefficient method which uses the uneven spaced 
spline coefficients to find the beginning, the peak, and the 
area of the lumps of a ship for classification. This method 
is better than the Fourier Coefficient method. The studyeee 


These methods is presented here. 
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lL. INTROBUGCEIGH 


It is quite important to classify the ships according to 
their types. Classification can be done in a number of 
different ways. The simplest is the visual method that is 
prone to error. Other methods for classification are the 
Fourier Coefficient method and the B-spline Coefficient 
method. In these methods the necessary information for clas- 
sification can be obtained from the superstructure profile. 

The Fourier Coefficient method samples the superstruc- 
ture profile at every chosen points. The function values at 
the sampling points are transformed into the spatial compo- 
nents. The logarithm of the magnitude of these components 
1s plotted and compared with the standard plot to recognize 
the type of the ship. 

In the B-spline Coefficient method, the spline coeffi- 
cients along the X axis and the Y axis are used to recon- 
struct the superstructure profile. The shape of the curve of 
the spline coefficients is in some way similar to the shape 
(the position of the lumps) of the superstructure profile. 
The ship classification may be achieved by recognizing the 


beginning, the peak, and the area of the lumps. 


sl 


mie ek See oes. ING 


Pee praocessitig~ss the procedure to obtain the superstruc- 
Eure profile of a ship from the IR image. Then, Fourier 
Meeerrcient Or spline coefficient methods can be used. The 


details of the preprocessing procedures are a follows. 


eee oaltA COLLECTION 


The data consists of the IR image of eight different 
types of ships. 
1. DD - Destroyer; "HALL" class. 


Zz. Container; The class is unidentified. 

co Freighter; The class is unidentified. 

4. AOR - Replenishment oiler; "WICHITA" class. 

5. LST - Tank landing ships; "NEWPORT" class. 

6. FF - Frigate; "GARCIA" class. 

7. CGN - Guided missile cruiser (Nuclear propulsion); 
"BAINBRIDGE" class. 

8. DDG - Guided missile destroyer; "CHARLES F. ADAMS" 
class. 


These images are taken from an aircraft which is flown 
at a height of 500 feet with a speed of 400 knots toward the 
Side of the ship. The aspect angles for these images are JO 
degrees which may be slightly off in some images. The inac- 
curacy of the aspect angle arises from the fact that the 
photos are taken while the aeroplane keeps on moving. All 
data of the images are stored in a digital magnetic tape 
with 256 by 64 bytes per image. Thus, the number of bytes 
required for each image is 16384 and each byte represents 
the intensity of a pixel. For each record of the image, a 
label is coded in the last 8 bytes as follows: 

1. Byte (16377) = Run number in each flight which passes 


i 


the ship. 

2. Byte (16378) = Video tape time code when the data is 
taken; in minutes. 

3. Byte (16381) = Video tape time code; in secomea 
Byte (16382) = Video tape time code; in thirtieth of 
a second. 

5S. Byte (16379) = Range in kilo-feet which is the 
distance measured fe snl the radar it may have an 
error 1 (Got sale ecu. 

6. Byte (16380) = Aspect angle; degrees from the bow 
of the ship. 

7. Byte (iSses) 
Byte (16384) 
testing. 


Ship class. 
ID, 1 = for training, Z, 3, 4, S325 


The run number and the time code together uniquely 
define a specific image that represents a single TV frame 
with no averaging. in. -aaedi ei] m, the time interval between 
the endof one image to the beginning of the other is 
approximately 1.5 seconds. Also, there are inherent random 
noise in the record which arises from the photo instrument 
and the process of storing them on to the digital magnetic 


Cape. 


B. SOBEE- OClepRalon 


The Sobel Operator technique is used to find the edge. 

To determine the edge, the Sobel Operator uses the differ- 

ence of gray levels of the pixels ina3 by 3 matrusee 
Shown in bigqure 2.12 

a, bee, a, 6, fnon and i are the values of the gray 

levels at the position Gf (x=l>y), (0 ot — 1. (xX, Woe 

(xtl,y), (-1,ytl1), (*%,yt1), and (2417742) respective ure 


Laplacian estimate is defined as 


Ot x Flxyy) ~ H+}, y) -[Fx+ Ly) Fx 29) 
Ox? 


GQ: 


14 


Rice DasisomYVeGceor for alle sdirections are (a-Z2btc), 
Meee), (a-2a4q), and (c=2f+t1)- Furthermore, each of the 
basis vector is convolved with the image as follows: 


along the xX axis 


dy = [F(x- LY=1) + 2flxy-1)+f(x+Ly-I)] 


Lee 
-[f(x-1, Y+|)+2F(x, y+) +F(x+1,y+1)| : 
along the y axis 
dy = [x+y +2f(+ LY) +F(x+1, ¥+ 1)] 
(ze ) 


~ [Hx Ly N+ 2t-Ly) +4(x-Ly+1) 


Since the magnitude of the resulting vector is the abso- 
lute value of the convolved results, the edge magnitude 
euey) ([Ref. 1], 


5) = (42442)! : (24) 


Note that the Sobel operator does not use the gray level 
at the position (x,y). The advantage of using a Sobel oper- 
ator over others is that the resulting edge is smoother due 
meoecms by 3 matrix approach. If we compare the Sobel oper- 
ator with the Laplacian operator, it is seen that the Sobel 
operator using the four basis vector as shown above will 
provide more accurate reading because of noise reduction in 
the original image. Hence, The Sobel operator is often used 


in the preprocessing operation. 


is) 





Figure 2-1 Sobel Operator. 


C. THE USE OF 4 Sebebeol kee 


When a Sobel operator is used at the edge of the image 
frame, the pixel level which lies out of the frame will be 
set equal to that of the adjacent pixel within the frame. 
The original image of a guided missile cruiser is shown in 


Figure 2.2. Since the result of the Sobel operator is numer- 


ically greater than 8 bits range, we have to rescale the 
result back to 8 bits range. This is achieved by deter- 
mining the maximum and the minimum of the gray level. They 


are then used to rescale the gray level in the results of 


the Sobel operator as shown in Figure 2.3. In some 
instances, the original image is very poor as shown in 
Figure 2.4. Attempts to determine the edge of this image 


failed as shown in Figure 2.5. 





Figure rZ az A CGN at a Range of 45000 feet. 
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Image from a Sobel Operator. 


a3 
Blured Image of a LST at a Range 67000 feet. 


igure 


a 





Figure 2.4 
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Image from a Sobel Operator. 
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Bagure 2.5 


DD. EDGE THRESHOLE = slo Clie 


Edge threshold strategies are used to extract the edge 
profiles from the Sobel results. "in this case 720s foe, 
the pronounced value of an edge element at x if g(x) is 


greater than certain threshold value [Ref. 2]. 


=g(xy) it g(x y)>> threshold 


G(x) 


(2.5) 
=O otherwise 


To increase the contrast of the image to’) a sii eee 


form, G(x7y) Isedertineadmas 


=255 if 9(x,y) Sthreshold 


G(x,y) 


. (2.6) 
=—@ otherwise | 


The choice of the threshold value is based upon the 
histogram of the edge image as shown in Figure 2.6. The 
estimated critical gray level is chosen so that a majority 
number of the pixels with value between 0 to 255 will fall 
below the critical value. Alternatively, histogram equa za. 
tion may be used to determine the desired threshold level as 
Shown in Figure 2.7. In this case, trial and error method 
was used in conjunction with the above method to obtain the 


threshold so that the correct profile is ascertained. 


1S. 





Figure 2.6 HiUSrocmanmwor Figure 2.3. 
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Figure 2.7 Cumulative Distribution of the Histogram 
in, Pioguce 2. or 
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PeeeeoW TO BATRACT THE PROFILE 


The edge image of the guided missile cruiser shows 
Miele Variations in the gray level as shown in Figure 2.3. 
These variations are caused by the noise in the original 
image. In this case, the choice of the threshold value is 
based upon both the histogram and the cumulative distribu- 
tion of the edge image so that it contains 90 % of the 
pixels. Therefor, the chosen gray level is 110 and the 


result is shown in Figure 2.8. 
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Bagqure 2.8 Siinonet GouOr sasecGN sin Figure 2.3. 


Pee Ow TO OBTAIN THE SUPERSTRUCTURE 


The original image. of the ship is taken from the aero- 
plane with different displacement from waterline to the 
peMperstructure in a rough sea, so that the profile of the 
ship with respect to the sea surface varies. Thus, we have 


Pomeliminate some information in the image of the ship by 


considering the superstructure only. In considering the 
everall ship structure, it is obvious that the largest 
distance is between the bow and stern span. Therefore, the 


bow and stern points are located first in the program as 
shown in Appendix A. Then we consider the slope of the bow 


and stern of the ship and set all the gray values below the 


Zak 


slope line to zero. Hence, it results in a superstruccureme. 
the ship as shown in Figure 2.9. However, in some images, 
there is a lot of noise in the background which cause diffi- 
culty in locating the bow and stern. Under these circume- 
stances, the dimensions of a ship are estimated by trial and 
error method. Then the gray level which lies outside is set 
to zero and an estimate of the bow and stern slope can be 


made. 





Puguinre we sUperstructure Profi levor ye maemo 


G. CONTOUR FOLLOWING 


The noise in the background of the original image, 
Yields wide edge structure. In considering this factoryiede 
contour following procedure tracking the inner part of the 
imagé in Figure 2.9 gives the Weuperserucemc o.oc ee The 
objective of this contour following is to describe =nemaaee 
and stern points of the ship, the direction, and the paces 
tion of the edge of the superstructure. The contour tracuag 
is done in the counter-clockwise direction which compares 


pixel value of O or 255 in a 3 by 3 matrix in the follows 


Nannie t:. 


22 


wi oMieeaoimeciace LErEMostmoolnt in the superstructure 
Meee In Figure 2.9 the contour profile tracing is accom- 
plished by examining the neighbors of a 3 by 3 kernel 
located at the curser position as shown in Figure 2.10. The 
SUrser is moved along the profile. All successive positions 
@meene curser constitute the contour profile of the ship. 
The testing procedure is explained below. 

ie Initialize the curser position to the beginning of 
the thresholded image with the gray level of 255 and 

the estimate direction of the next position. 
feeeeeneck for reaching the end of the profile, if it is 
Hee ele emailing) Of Enews prorile then stop, if not go 


iso: 3. 
oe Check for the estimate to see whether it is in the 
Smeeectemon of, North. Bast, South, or West. Tf the 


Seebecetom tt 16. Norehn then go to 4, if, 1 CN aS east 
iC lovey ele) key fs) - Pim omootenmchen Go to 6, and 1f it 
is West then go to 7. 

4. Determine the next position with the gray value of 
255 in the counter-clockwise direction as shown in 
Figure 2.11. Store the position found and move the 
eurser -€6 that position. Update the estimate 
SerheOecrTom EO Ene one Vase £ound- then go to Z. 

epee ihne procedure is the Same as that in step 4 except 
search pattern is shown in Figure 2.12. 

6. The procedure is the same as that in step 4 except 
search pattern 1s shown in Figure 2.13. 

7. The procedure is the same as that in step 4 except 
search pattern is shown in Figure 2.14. 

The flow chart of the procedure is shown in Figure 
eo, and the detail of each procedure are included in 
Appendix B. The testing procedure have to be performed in 
such a maner that the resulting contour is a good represen- 


tation of the superstructure line. The result of the contour 


Zo 


image 15s shown in 220gure Zoe If the superstructure in 
Figure 2.9 has wide edge, we can not find the contour image. 
We have to use the additional step which is called the 


"Closing" operation. 


Figure 2 sk® The 3 by 3 Kernel. 





Figure 2.11 The Step Checks in the North Direction. 
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Peguine 2.12 The Step Checks in the East Direction. 





Regure 2.13 The Step Checks in the South Direction. 





Figure 2.14 The Step Checks in the West Direction. 
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Figure ZAG Contour Image of a CGN at 
a Range of 45000 feet. 


Pec LOSING OPERATION 


Some results from the superstructure extraction process 
are discontinuous because the gray level of those areas of 
the structure is less than the threshold value. If we 
decrease the threshold value, the details of the superstruc- 
ture are effected. It is necessary to use the "Closing" 
operation which consists of the "dilation" process to smooth 
the superstructure profile used the "Erosion" process. All 
direction are dilated similarly which causes an smoothing 
effect on the edges, The superstructure increases in total 
area. Then, use the "Erosion" process to shrink (subtract) 
the dilated part in all directions, thus obtain the smoothed 
eeeerstructure with the appropriate size. ines, Closing 
Beeocess employs the dilation process which yields an output 
image from an input image. The Dilation results is shown in 


meoure 2.17. 


Ze 


O orinal image at gray level=255 
X center of kernal original 





FLOUrEG 2.7 The Process of Dilation. 


We examine the gray value of each pixel in the original 


image. Begining from the pixel at the first row and the 

first column. The procedures are the following. 
1. Considering one pixel in the original image with the 
kernel (B) centered there. If at least one pixel in 


the kernel has a value of 255, we let the gray level 
of the output image at the center of the kernel be 
Za 

Zz. We shift the center of the kernel 1 column to the 
right. Then following the same procedure as in step 1 
until the last column is reached. 

3. We shift the position of the kernel to the next row 
and starting from the first column. Then, followan 
the same procedure as in step 1 and step 2 until the 
last row and the last column is reached. 

The result obtained from the dilation process is “an 
image with enlarged structure. The second procedure in the 
Closing operation is the Erosion process. The Erosion 
process perform the same procedures as the dilation process 
GxGepu fer oreo =. If every pixel in the kernel are 255, we 
let the gray level in the new image at the center of the 
kernel be 255. Otherwise, 1t will be QO. The output image 


obtained will have smooth edge with minor change occurring 


Zo 


Gm the edge detail as shown in Figure 2.18. In this case, we 
use an kernel of size 3 by 3. If we increase the size of the 


kernel, the details of the image are decreased. 





Eeogure 2.18 The Profile after Dilation and Erosion 
(Closing Process). 


Peer OrlLE ROTATION 


Often in the contour image, the first and the last point 
of the superstructure are not at the same horizontal level. 
We have to rotate the contour image by setting the two 
points to the same horizontal level. HPOwereOenevate 2c from 
the Y, X axis to the Y', X' axis is shown in Figure 2.19. 


If the angle value 9 1s positive, it will be 


xX! _ cos -sing x 
baie <= sing cos Be 
If the angle value 9 is negative, it will be 
| cos sin Xx 
a -sin§ cos vy 
Bereesome OF che Contour profile, part of the profile 
after rotation will be out of the image frame. Then we have 
Bemeeshiftt this contour down by 20 pixels position. The 


rotated profile is shown in Figure 2.20. 


Jd 





Figure 2.19 Rotation Emocessc. 





Figure 2.20 Rotated Prof 
a Range of 45000 
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Poe woOURTER  CObmare 1 ENT Mz THOD 


Wee have obtained the ship profile in the previous 


chapter. ReGwewrmace feattses OUL Of this profile for clas- 
sification purposes, we will use the Fourier Transform 
method. 


fie FOUrIGr  Cransform of Che ship profile showed that 


the transform coefficients depend upon the ship's dimen- 


sions, its superstructure, and the distance between the 
Camera and the ship. ceo meonorlmhcer(s.) lS) va discrete 
Mam@eclOn with 128 sample points. Weve Cols iieisis Jd @ lis allie ie 


jransform can be written as 


M-/ 


(x)e x p(-j27ru x/N) (3.1) 


meme Ww, KX = O, Ll, 2, «ss... 


samples. 


, N-l. N is the total number of 


Tf the direct calculation of the discrete Fourier trans- 
Mmeemeets Chosen, the number of complex multiplication and 
Meewetton will be equal to N; i.e. to obtain F(O) would 
require complex multiplication and addition N times. In 
order to reduce the computation, we use the fast Fourier 
Mmemamstorm algorithm. Thule ew weduatki1on 3.1 (Ref. 2] can be 


separated into even (u) and F ogg (u) 


nia U xX 
Feven(U) “Lh (2x) ae 
Moe 


(342) 
at U x 
Fadl) = 1) f(2x+1) Wi, (3.3) 
X=0 


Sa 


W =exp(-j27t/M) (3.4) 


M=N (335) 
Z 
U 
rtu)= 2 Fev ent) + gq (tt) Wow | (3.168 
a U 
F(u+M) = [Feve nll) Fog aU) My] (3.7) 
Using this method, we have reduced the total number of 


complex multiplication and addition to Nlog N. In thisWeaars 
we have N = 128, thus, the total number of complex multipli- 
cation and addition w2tl bese sc. 

First, we divide the rotated profile image into 128 
division se Since the distance of each division is equal to 
the amount of the pixel between the bow and the stern of the 
ship divided by 128, we use the distance perpendicular from 
the horizontal line between bow and stern to the highest 
point of the superstructure as the sampled values. The 
result of the fast Fourier transform is a complex number. 


Then we use 


M(u) = 10 g[ 1+ G(u)] (388 


G(u) 1s the magnitude of F(u). A value of 1 is add to the 
magnitude to avoid negative logarithm result, the results 
obtained are as shown in Table I and 


i Figure 3.1 - Destroyer 


Z. Figure 3.2 - Container 

oop Figure 3.3 - Freighter 

4. Figure 3.4 - Replenishment oiler 
ae Figure 3.5 = Tank landing esas. 


Sz 


S.. Figure 3.6 - Frigate 

Yeme Figure 3.7 = Guided missile cruiser 

8. Figure 3.8 - Guided missile destroyer 

Mm the results minor difference of the shape of the 
Fourier coefficients can be noticed visually. Beet. as 
difficult to implement a program to detect these minor 
difference in shape. Therefore, a Second method was used to 


handle this problem. 
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Logarithmic Magnitude of the Fourier Transform of a Freighter. — 


Figure 3738 
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Logarithmic Magnitude of the Fourier Transform of a AOR. 


Figure 3.4 
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Logarithmic Magnitude of the Fourier Drains BO sof amino ls 


Figures 35 
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Logarithmic Magnitude of the Fourier Transform of a FF. 


Figure 3.6 





40 


Logarithmic Magnitude Of thestotmtcr Tranctorm of waeon. 
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Logarithmic Magnitude of the Fourier Transform of a DDG. 
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Destroyer at range 77000 feet. 


DD 


Cointainer at range 28000 feet. 


Cont 


Frig= Freighter at range 40000 feet. 


Replenishment oiler at range 78000 feet. 


= 
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AOR 


Tank landing ship at range 51000 feet. 


— 
= 


Loa 


Frigate at range 49000 feet. 


FF = 


Guided missile cruiser at range 45000 feet. 


GG 
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Guided missile drestroyer at range 41000 feet 


—= 
= 


DDG 
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Peeve i ATION OF SHIP SUPERSTRUCTURE WITH RANGE 


Same of the practical problem of using the ship profile 
Memenat iC 1S Sensitive to range variations. Close ship 
profile has more details than the far away ship profile. 
The dependency of the geometric size of the profile with the 
range is discussed in this section. 

Assume that the object is centered on the camera axis. 
The field of view of the camera, the number of the pixel in 
the image, the size of the image, and the size of the object 
are known. Our problem is to determine the distance between 


the camera and the object. 


a. Background 


The camera system is similar to a human eyes system. 
The light reflected from the object goes into the eyes. The 
image of the object falls on the retina, the signal is sent 
to the brain in some electrical from, and the brain changes 
it to a from that human can perceive. But, in the camera 
system film or senser are used to pick up the image. The 


function of a camera 1s shown in Figure 3.9 





Pacuires. 3.9 The System of the Camera. 
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Fi GQure 76.510 Range Calculation. 


‘The distance between the lens andthe image plane 


can be adjusted in order to have a clear image on the film. 


The camera has a field of view (FOV) angle as shown in 
PIQUE Bes] The object has to be in the field of view of 
the camera. The determination of the distance is shown in 


Figure 3.10 
For simplicity of calculation the inside cilwien- 


camera is flipped to the same side as the object as shown in 


Figure 3.10. when the angle of the FOV isCY, I/2 is the 
half of the full image size, D/2 is the half of the dimen- 
Sion of the object, XK is the distance from the lens tomeme 


image, and R is the distance from the lens to the object. 


Assume that the’ distances 2 6 distaneommrm 2e and 


angle(¥/2 are known. Then, the distance R can be determined 
by 
MS = Apa (al ave? 
rans ag 
tand =] 
2 ax ee 
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an (ly) 
ZL Sale 
2 


Assume that the angle resolution of the pixel of the 
field of view of the camera is equal to 0.2E-3 radian per 
pixel. The number of pixel of the frame is equal to 256. 
The size of an image is I pixels. The dimension of the 


object D in feet is known. The field of view in angle is 


OF SG) Zle= Sy 2 Se 1D) 
X=2560 qd 
D tang Coe.) 


Mmfere d is in unit of pixel. 





tan@ =1 d 

2 2% (3.14) 
ton = 1 tan 

a LAG 2 (oer Ss) 
R=D | 

2 Tancy' Eels 

Z 

R=128D 

ial nO (ODSTE (3.17) 


When the length D of the object is known, from equa- 
mom 5S.1/ we Can estimate the distance from lens to the 


object and is shown in Table II 
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TABLE 


Range Estimation 


| | | yy | | ; ml 
| CLASS | I(pixel)] D (ft) | R (k£t) | R(kf£t) | (R=-R)100] 
! | | | | | -=aiaie | 
|MEASURE | KNOWN  |CALCULATE|RADAR DIS] R | 
| | | | | | | 
| - anno nn nn nn nn nn nn nn nnn nn enn nnn nen nnn nee -2---- | 
| | | | | | | 
| DD1 | 96 | 418 [Pe eG el ea fy releesr ees 
| | | | | | | 
| DD2 | 80 | 418 | 26.119 | 85 Wess: 2 el 
| | | | | | | 
| AOR] © |, “107 "ese sei Wc A | 60.53 | 
| | | | | | | 
ib RORZ.) SI 96 | 659 (340505) see | 59.63 | 
| | | | | | | 
| LST. 9] 176 "We 2203 eles een ees i) Soviets, 
| | | | | | | 
| LSt2. ||) “134 9) @522.535 ioe oie | 65.82 | 
| | | | | | | 
| CCN] | 147 #+2'| 565 iol isos | ars (357-3 a 
| | | | | | | 
lC6h2- (| 6119 Seieesos | 23.934) ese | §6.85 | 
| | | | | | | 
\ BBG] |. Wee mess 17 937 - a eee | 63.11 | 
| | | | | | | 
| DDG2 | 90 | 437 245 24 7c 2 | 62.08 | 
| | | | | | | 


DD1,DD2 = Destroyer at range 79000 and 83000 feet. 


AOR1,AOR2 Replenishment oller at range 78000 and 88000 feet. 


LST luke 


Tank landing ship at range 51000 and 62000 feet. 
CGN1]1,CGN2 = Guided missile cruiser at range 45000 and 64000 feet. 


PDGT, DEGZ 


Guided missile destroyer at range 41000 and 64000 feet. 


The error distance between the estimated distamee 


and calculated distance in percentage is ((R - R)/R) 1008 


Zz Remark 


Calculated distance error in R may come either from 
the pixel measurement in an image or from the angular reso- 
lution estimation. The distance that is stored in the image 


label has an error of about 1 to 2 kilo-feet. If the angle 
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of the field of view is accurate, we can estimate the number 
of pixel in an object image. Then we can classify the object 
by comparing the number of pixel of the original image with 
that of the test image. Some known system parameters can 
help to determine the range of the target ship. The problem 
is that existing errors in the system parameter causes in 
the larger errors in the estimated range. Consequently, they 


are not very useful in classifying the ships. 
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IV. B=SPLINE COEFFICIENT ET HOD 


Using B-spline coefficient is another method to describe 
a ship profile. This method uses the B-spline coeftficivemes 
EG determine the beginning, peak, and area of lumps which 
contain significant information about the type of Che jshujoes 
The comparisons of the knot position (in parametric value) 
from the midships to the peak or beginning of the lump, can 
be very helpful. Different ships will have different lump 


positions and sizes: 


A. BACKGROUND 


A spline function 1S a piecewise polynomial used to 
interpolate points. This kind of curve is smooth and the 
discontinuities in its k-th derivative is as small as 
possible. In this case, Cubic spline was used, where the 
first and second derivatives for any set of interpolating 
points are continuous, while the third derivative may be 
discontinuous. The reason for using Cubic spline is to Keep 
the third derivative discontinuity as small as possible and 
the curve as smooth as possible. The B-spline calculation 
procedure is also very stable. In our case the order of 
spline used is 3, while the 1-st and 2-nd derivative are 
continuous. The choice of 4-th order spline function is due 
to the fact that it is generally sufficient for most sie 
profile curves. Discontinuity, in a sense, may be stated as 
the jumps of the third derivative, which is the means to 


control smoothness of the connecting pieces. 
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Ere o-oPLINE APPROXIMATION WITH FREE KNOT 


As mention earlier, the B-spline function is used in the 
spline approximation with free knot. The free knots is some- 
times called uneven or lrregular knots; Ghaen1ec - mo fixed 
number of knots are used. Furthermore, the spline position 
need not be on the original curve so as to minimize the 
number of knot position while preserving most of the details 
@mecnie Original curve. 

Perst, Minimize the value of the lack of the smoothness 
Weeyecefined as [Ref. 3] 


| 
anh Cals 


Mme C1 is a coefficient of spline at the Knot position. 


aij is defined as follows 


= M Ady +0) - Mi, Met 1) 


(4.2) 
where Mi kel is the normalize B-spline function and defined 
) 
as 
M 2h, Gust tx 
spend ed) a i+ k+/] p t (Heyes i ney 1 9KM ) 
(4.3) 
and G,(t;x), t are defined as follows 
k k : 
=(t-%), = (f=) if t2>x 
g,(4;x) (4.4) 


=O eee 
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Ds ap Aes pene Gia DECC) stands fou Che l-th® divided 
difference of the function f(t) on the point 4). 
where t is the position values of knots in term of Z param- 


eter. Z parameter is defined as follows 


2 2,2 
Z(L) = 21-1) + [(x(1) = X(1=1)) + (¥() - Y= 0) 75 (4.5) 


Z(O) =O 
(4.6) 


where I is the number of the sampling point Wena 
S12 Seale 


Second, the smoothing is subjected to a constraint 


O(c) < s (4.7) 


where S is the smoothing factor, O(c) is the weighted sum of 


the square residuals defined as 


tc) = Suly-F 
io = Dil ~ CUM 44/9) (4.8) 


Ay, Vy are the values of X and Y at Z parameter of the 


sampling points; We 1s a weighting factor for all sampling 


“points [Ref. 3] defined as 


-2 
Wi = (Oy,) (4.9) 


where Wi 1S an estimate of the standard deviation of Yj 
Then, the value of S is in the range m+V2m. If nothingmame 


known about the statistical standara wi eevt ae onver ty each 


310) 


W can be equal to one and the S can be determined by the 


metal and error method. 


ab. ove sve ore Leealieio, se vere! epson Mak .Om. USiing B-spline 
le bheke = lies! 


There is a difference between the interpolation and 
the approximation. In toeewoolacion;, “the number of Knot 
positions required are equal to the number of sampling 
MPemmes and the value of S in eq(4.7) is small. In this case, 
the computation time required will increase tremendously. 
feeweras, in approximation, it is not of great concern that 
the approximated value at a position be the same value of 
the sampling point, and most of the information still remain 
in the original curve. Thus, decreasing the value of S will 
result in the large number of the knot positions and the 
final curve obtained will be similar to the original curve. 

The justification for using approximation rather 
mga interpolation approach is that though the resulting 
curve may not be the best fitting curve, 1t is smooth and 
close enough to the original. The approximation spline func- 
tion has less number of knots than the number of samples, 
which reduces the total processing time. In approximation 
approach, when the appropriate choice of S value is made, it 
fae in turn, Semeractersurticicnemmumoer ot Knot positions 
required to provide a close approximation to the original 
‘curve. The ratio of sampling point to spline coefficient is 
10 to 1 as shown in Figure 4.1. ia picure 4,1], a guided 
missile cruiser ship at a distance of 45000 feet, with 290 
Sampling points and the associated B-spline knots are shown. 
The original samples are plotted in solid curve. The Knot 
position are show by small circles in Figure 4.1. After 
B-spline approximation, the number of the Knots are reduced 
to 33. Hence, this technique is essentially a kind of data 


compression technique. 


orl 





Figure 4.1 Plot of the Original Data and the Approxaieee 
with Free Knot. 


C. TO DETERMANE THE * KNOT POST lethiyy At THE B=-SE ees 
COEPETCIENAS 


The approximate smoothing factor in eq(4.7) is to be 
calculated before using the subroutine "PARAM" [Appendix C]- 
In the subroutine, there is a check on S factor every time 
it is run. If the S input values do not satisfy or mecumene 
criteria, the program will vrecurn with a error code IER. 

There 1S a parameter NEST which is set to a constant 
value. This value determine the dimension of the knot posi- 
tions which relates to the array T(NEST), Cx(1..NEST) me 
Cy(l..NEST). The NEST is an overestimate of the dimension of 
the arrays set by the user. The limitation on the value of 
NEST for subroutine "PARAM" are as follows [Appendix C] 

is. 2k NESs Make) 

2. Typically, value of NEST M/2 
where M is the number of the total sampling points and k=3 


1s the highest order in the B-spline function. 


The subroutine "PARAM" produces several outputs as, N 
Che total number of Knot pesweroms- T(N) the array of the 
value of knot positions in Z parameter, Cx(N) and Cy(N), 


B-spline coefficient of X functions and Y functions for ime 


positions defined in array T(N). 


SIZ 


ime Limitation On B-spline Coefficient Determination 


ti theme rvaliecwor NEG@ecmnetoorsmall, the user will 
receive error code IER and the number of knot positions will 
appear to be dense in the first part of the curve. While 
Sparse in the other part. The example of an incorrect selec- 


fmeem of NEST is shown in Figure 4.2. 





Peeare 4.2 Knot Selection if the NEST Parameter is too Small. 


Another problem which causes difficulty in program- 
ming is that the main program is in PASCAL while the subrou- 
tine is in FORTRAN. As for the PARAM program in FORTRAN, the 
array index value starts from 1, while , for the user PASCAL 
program it starts from 0. Therefore, in linking the main 
program to the subroutine, one has to keep in mind of the 
difference. in, acd t ton, the FORTRAN programming logic 
Seructure 1S so complicated that, when an error occurs, it 
Meeevery difficult to debug and locate the error. 

The values of Cx and Cy depend upon uneven knot 
Mertetons, and they contribute controlling effect to the 
mecoOnstructed curve which would be both smooth and close to 
mie Original curve. In running the B-spline approximation 
program for the first time, the values of Cx and Cy of the 
last 4 knots at the right end are zeros. However, in running 
it again, PiGrweacse thesValuc Of.5 did not yield zero values 
feeee wana Cy at those points. This, nevertheless, has no 


PmEect On the reconstruction of the curve. 


Sy 


Plots of the reconstructed profile of tnree shipewee 
the CGN class using different value of “SS, result in 
different number of knot positions and are shown for compar- 
ise in Figure: 4.37 4.4 and Figure 4.5. The original 
Samples are plotted in solid curve, the reconstructed curve 


is plotted in dashed line in Figure 4.3. 
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Figure 4.4 A CGN at Range of 55000 feet. 
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Figure 4.5 A CGN at Range of 64000 feet. 
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2. Selection of Smoothing Haceor Valuci ey 


aw 


It is found by trials and errors using the compijges 


program PARAM that in order to retain maximum information of 


the profile curve, the number of knot positions required 
should be in the range. of 25 Commo. The number of Knot 
positions depends upon the value of S which must be set 
accoOneaingly. The appropriate choice of S$ to satisfy the 
condition stated above, is important. For the class etme 
guided missile cruiser(CGN), three ship images at 3 
different ranges were selected. Then, run the appropriate 


program for various S factors to see how the number of 
knot(N) will vary. The results are shown in Figure 4.6. 
Plots of N vs S in Figure 4.7 through Figure 4.13 show Gia 
in most cases, the value of N decrease quite rapidly when 
the value of S 1S an the range cl Oeaercer! ce, and gradually 
for S factor in the tance of 008-2 oo, thereafter, the 
value of N decreases very little. Obviously, the curve seems 
to decay exponentially. Furthermore, for some classes of 
ships changes are more pronounced than the others which is 
probably due to the actual number of Knots present in the 
profile. For guided missile cruiser ship(CGN) with Z2 Puimeee 
the number of knot positions required can be 33 as shown in 
Figure 4.6 We select the factor S to be about 100. 

The selection of the factor S depends upon the 
number of the original sampling points. If the factor Sie 
small, large number of knots are needed. When the factor S 
is large, small number of knots are needed. When the number 
of knot and the Cx and Cy coefficients are small, the 
B-spline coefficients Cx and Cy obtained can not be used to 


reconstruct the curve close to the original profile. 
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Bree tbe OulcepuUE @x,Cy and Knots Profiles 


From the previous study the factor S is selected for 

each of the class of ship as follows 

1. Destroyer(DD), S = 20.0 
Container, S = 100.0 
100.0 
Replenishment oiler, S = 20.0 
Meonkwwuanding ship tLst), Ss = 25.0 
Frigate(FF), S = 100.0 


Guided missile cruiser(CGN), S 


Freighter, S$ 


172 SO 
Guided missile destroyer(DDG), S = 125.0 
The plot of the B-spline coefficients,Cx and Cy,and 


On vn WM SP W ND 


X,Y at the positions of the sampling points vs the Z@ param- 
eter are shown 1 Jol Figure 4.14 through Ba ule e472. ie 
Observation and comparisons of the curves show that the 
values of Cx and Cy exhibit changes similar to that of X and 
Y except that the variation of values leads that of the X 
aad Y. This is due to the fact that Cx and Cy have to act 
Memmeene Controlling factor for the reconstructed B-spline 
curve to get the result close to the original curve. 


Examination of the plots of the results of X and Y 


show that when X is increasing monotonically, Y is almost 
@emstant; but when X is almost constant, Y is increases or 
decreases. This behavior Ceohabl VG tO MrEorlle reconstruction 
may be explained as follows. For a ship profile when X is 
increasing and Y not increase too much, this may be inter- 
preted as an almost leveled profile. When X is) almost 
constant, and Y may be increasing or decreasing, it may be 


interpreted as the beginning or the ending of the lump. 
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Figure 4.14 


and Cy vs Z for a FEF at a Range of 43 K-ft. 
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Plot X,Y,Cx, and Cy vs Z for a CGN at a Range of 45 K-ft. 


Figure 4.16 
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Figure 4.17 


and Cy vs Z for a DDG at a Range of 41 K-ft. 
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and Cy vs Z@ for a AOR at a Range of 78 K-ft. 
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Figure 4.21 


at a Range of 28 K-ft. 


BP 


ELons 


SHLP DESCR Varro: 


The shape of the ships depend upon the shapes and posi- 


of the lumps. Characteristics of the lumps Siren 


different type of ships are as follows: 


ty: 


Frigate - The beginning of the lump is at 1/6 of the 
ship length to the left side of the midships and 
the peak of the lump is at the mast which is located 
at the midships. The termination of the lump 
is approximate 1/3 of the ship length from the 
midships to right side. In AaAdeai tion the average 
height of the lump is a little higher than the 
level between the bow and the stern, and its size 


is 1/2 of the ship length as shown in Figure 47ze 


“go os ie < ® ww > ne) 

Pi Damn, SOS i wet in tar 94 : 4 - = . 
oa wee TL EALAGK ce 

PRES Fr FTO te bese ae by oo? od Oh) Wie 
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Figure 4.22 Frigate. 


Container - The beginning of the lump is at 17am 
the ship length to the right side of the midsimiee 
while the lump is high and terminate at the stern. 
The lump appears to be in a rectangular shape with 
small crane. 

Tank landing ship(LST) - The beginning of the lump is 
at 1/4 of the ship length from the midships to the 
left side; the height of the lump is higher thaneenae 
level between the bow and the stern by a small 
margin, while its highest point is located at 


approximate 1/6 of the ship length from the midships 
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EOuMGnoulepe wc ae wither small difference from the 


average Nigh as shown am Figure 4.23. 


NEWPORT, “NEWPORT” Class LST USA 120} 





Figure 4.23 amie sania ng sol itp( EST): 


Guided missile destroyer - the beginning of the lump 
starts at approximately 1/4 of the ship length from 
Ehe Mlaships toeene Left wside with its highest point 
at the mast. After this the slope will decline ina 


very rapid fashion fee kowAl MiG) a noticable deck 
distance, then the beginning of the second lump 
occurs due to the redome presence. Therefore, the 
lump will be narrow with great height, speveud aia) 0a 


terminate at approximatly 1/6 of the ship length from 
the midships to the right side. Furthermore, there 
is a Latte lump near the stern, this 
distinguish destroyer of the same size as shown in 
Figure 4.24. 

Peceeover-luMp Characteristic will be similar to that 
of guided missile destroyer except for the small lump 
as in Figure 4.25. 

Guided missile cruiser - The beginning of the first 
im 1S amy iz or the ship tengen from the midships 
to the left side with highest point at the mast. The 
lump size is large both in length and height and its 
height decreases to the point which is a little above 
the level between the bow and the stern. Therefore, 


the second lump begins with almost the same size as 


Ta 


"CHARLES F ADAMS: Clace NOG Aust (3), USA (23) 
f PERKIN Class sinvuiar with deca nose between funnets)} 





Figure 4.24 Guided Missile Destroyer(DDG). 
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Figure 4.25 Destroyer. 


the first one. The second lump ends at approximately 
1/4 of the ship length from the midships to the gre 
Side. Furthermore, the distance between the peak of 
both lumps is less than that of the Guided missile 


destroyer shown .in Figure 4.26. 
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Figure 4.26 Guided Missile Cruiser(CGN). 


Wes 


7. Replenishment oiler (AOR) - There are 2 little lumps 
with the first one starting at 4/5 of the ship length 
from the midships to the left side exhibiting rapid 
decrease in height to the level between the bow and 
the stern. The second lump starts at 1/4 of the ship 
length from the midships to the right side with slow 


decline to the point near the stern. 


8. Freighter-with 3 noticeable lumps, the first two 
represent lifting crane with the first high lump 
approximately 1/6 of the ship length from the 


midships to the left side but narrow; meanwhile, the 
second lump is located at approximately the 
midships with the same size as the first one. The 


beginning of the third is at approximately Se © 
the ship Length L£rom ithesmicdships to the right 
Side with large size. Kemi S higherwrenan the first 
two and with gradual decline to the stern. 
From the lump characteristic of different type of ships, 
we could distinguish the type of a ship based on the rotated 


Superstructure. 


Bee ool P CLASSIFICATION 


PMmEicomcanc aspect In GCategorizing types of ships is to 
recognize the lumps above the main deck. Therefore, dee eis 
‘useful to plot the positions of X and Y vs the Z parameter 
Where the Z parameter can be determined from eq(4.5) where 
Z2(I) 1s an array (1..M) as shown in Figure 4.27 where X is 
ploted in solid line and Y is ploted in dash line. Then, 
plot the B-spline coefficients Cx and Cy vs T(N); the knots 
position in Z@ where N is the total number of knots as shown 
in Figure 4.28. The Cx is ploted in solid line. The Cy is 
ploted in dash line. The Values of Cy will be close to the 


curve of Y while the values of Cx for the same knot posi- 


TES, 


tions as Cy is exhibiting monotonic increastigeen wee When 
the difference in Cx is decreasing or zero, the value of Cy 
will have a pronounced change where the beginning or the 
ending points of a lump can be detected. The value of knot 
position(T) at those points related to the sizes of the 
Cx,Cy can be determined. Finally, with the values of Cx and 
Cy at those points known, the area of the lumps can also be 
determined. 

Hence, from the ship's characteristics and information 
derived from the above procedure, classification of ships 
can be made by considering 

First, the number of lumps detected, 1, 2, or 3. 

1. The number of lump=1: Frigate, Tank landing ship 
2. The number of lump=2: Destroyer, Guided missile 
cruiser, Guide missile Destroyer, and Replenish- 


ment oiler(AOR) | 


3. The number of lump=3: Freighter and Container 
Second, the position of a lump relative the midships is 
measured. This quantity is scaled by the total length of 
the shape This scaled quantity will be invariant with 


respect to the different ship sizes at different ranges 
He lausge|y the area of the lump is normalized to the ship 
length squared. 


i Lump Detection 


As shown in the plot of X, Y, Cx, and Cy vs Z@ param- 
eter in Figure 4.15, when the difference (ACx) between 
successive value of Cy varies from increasing to decreasing, 
and then, to increasing again, Cy exhibits noticable varia- 
Elon: From observation of Cx values, it is seen that the 
difference (ACx) always has variation in the same sequence 
as stated above. Therefore, only Cy values are taken into 


consideration in the program that detects lump. 
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There are two distinct procedures used in dealing 
with two different kinds of lumps big or small. Thus, it is 
necessary to distinguish in the first place whether the lump 
is big or small. For the big lump, the differences of Cy is 
meeomerve for all initial four knots. It continues to the 
peak and decreases toward the ending of the lump. For the 
small lump, the difference of Cy is positive for the first 
im OmRtOotcs and stay constant or positive for the third knot, 
Peuemeric difference of Cy for the forth knot is negative, 
thereafter, the program proceeds to establish the following 
values for each lump detected: 

First, the knot positions (Z value) at 

1. the beginning of the lump 

2. the ending of the lump 
Second, the Knot number for 

1. the beginning of the lump 

2. the ending of the lump 
Third,Number of lumps detected. 
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a. Procedure to Detect the Lump 


Check a big or small lump by testing the conditions 
of three varying values of Cy increments (ACy) in 
sequence. If they are positive it is a big lump and 
set the flag-lump to 1, otherwise it is a small lump 
and set the flag-lump to zero; then go to 2 

Check the present knot position to see where its 
Z value equal to zero or to maximum Z value, if it is 
4 Maximum then (Stop, 1 Gstice Go eos 

If the process begins at the first knot position, set 
the begin and the flag-end to zero; check the status 
of the flag-lump for 1 (big lump) or O (small aun 
if it is a big lump, then Ge coe. If it is a smaue 
lump, chen Ge =teose.- 

Check the status for the begining or ending of the 
lump. If the flag-begin is l, it represents that the 
beginning of a lump is found, then go to 5. Other ae. 
it S met. found) sehen vce —sromer 

Find the ending of the big lump by testing the 
conditions of 3 values of Cy increments in sequence. 
The first 2 should be negative and the third should 
be “Constance, of positive. If the condition are 
satisfied, store the Zvalue of the position of the 
third knot, and the flag-begin to zero; then Gomes 
LO: 

Find the beginning of the big lump by testing ier 
conditions of 3 different values of Cy increment 
(ACy) in sequence. The first Cy should be negative 
OF OnsStant ; the second should be positive, and the 
third should be positive. I£ the conditions ia 
Satisfied, store the Z@ value of the position of the 
second knot, and set flag-begin to 1; then go to ie 
Check the status for the beginning or ending of the 
Pumm. If the flag-begin is 1, it represents the 
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ae) 


beginning of a lump is found, then go to 8. Otherwise 
iin snot founcs chen "qouto o. 

Find the ending of the small lump by testing the 
Senglei1ons Otmommvalues or Cy increment (/\ Cy) in 
sequence. The first one should be negative or 
GOnstant.- If the conditions are satisfied, store the 
Z value of the position of the second knot, and set 
Ene wt lag=begqin toe Zero; then go to 10. 

Find the beginning of the small lump by testing the 
conditions of 3 values of Cy increment (A Cy) in 
sequence. The first should be constant or negative, 
the second should be positive, the third should be 
constant. If the conditions are satisfied, store the 
Besteienm Of the second knet, and set the flag-begin 
ipG? eamclien Gosto 10, 

Move to the next knot, then go to 2. 


This procedure is shown in Figure 4.29 andthe detail of 


each procedure is shown in Appendix D. 
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mm rr a rm me mem 


iie areasumader the lump can be determined by 


AAREA = Ads cy + CyACy (4.10) 


Suppose there is a lump as shown in Figure 4.30. 
The area increment can be calculated by using eq(4.10). Then 
the area under BC(which is negative) is added to the area 
under AB. 





Figure 4.30 First Procedure to Determine the Area. 


Next, area under CD is calculated, which is positive 
pmramacda this to the last resulting area. Obviously, the sum 
would represent the total area of the lump as shown in 


Figure 4.31. 
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Figure 4.31 Step by Step Procedure to Determine the Area. 


FEF. CONSTRUCTION Of Fee beet tener: 


The characteristics of ships used for classification 
are, the number of lumps, the area of lumps, the knot posi- 
tions(Z value) and where the beginning of the lump and lump 
maxima relative to the midships are located. 

In view of the above characteristic, 1t is necessary to 
find the relationships among them to make classification 
possible. Therefore, for each of the eight different class 
of ships, the decision tree 1s constructed which 1s based 
upon relationships observed in the plots of the knot posi- 
tion (2 values) for the beginning of the lumps normalized by 
the total ship length (Z value) VS. the number of lumps; the 
area of the lumps normalized by the total ship length 
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Mampvalue) squared vs the number of lumps; and the Z value of 
the lump maxima, as shown im. FPaiqure 4.32 egusewwleta 
Figure 4.39. 

Thus, according to the number of lumps presented in the 
Peerile used as the first criteria, ships can be divided 
moors GAistinct groups. Then, classification, can be made 
by comparing further characteristic as shown in Table III, 
and the complete decision tree constructed from Table [III is 


Shown in Figure 4.40. 
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Figure 4.33 Plot the Beginning, Area, and Peak of a LST. 
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Figure 4.34 Plot the Beginning, Area, and Peak of a DD. 
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Figure 4.35 Plot the Beginning, Area, and Peak of a DDG. 
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Figure 4.37 Plot the Beginning, Area, and Peak of a AOR. 


o5 


BEGIN FROM MIOSHIP 


= Oso 


AREA QF LUMP 


Oe) 


0.005 





NOe OF eiUee 


FREIGHTER AT RANGE 40000 FEET 
$=100, N=S34 


— 
it 


2 = FREIGHTER AT RANGE 47000 FEET 
§=50, N=26 
fo oe 


FREIGHTER AT RANGE 53000 FEET 
$=50 , N=24 





/ > 3 NO. OF LUMP 


AREA OF LUMP 


0.0! 





—0.5 0.5 FEAK FROM MIOSHIP 


Figure 4.38 Plot the Beginning, Area, and Peak of a Freighter. 
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G. SUMMARY 


The decision tree which 1S constructed from Tablemme 
does not provide l100-percent correct classification fon ga 
types of ship images. Furthermore, the presence of noise in 
images may cause complicatations in classifying the ships. 
For example, with excessive noise in the original image, 
Sobel operator for edge enhancement in the preprocessing 
process, Still yield result with residual noise present in 
Figure 4.41. These residual noise 1s undesirable since it 
causes failure to extract profiles which retain necessary 
informations from the original images. The subsequent clas- 
Sification of profile by the Fourier Coefficient method and 


the B-spline Coefficient method becomes difficult. 
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Figure 4.41 Noisy Image. 


As discussed before in order to reduce the noise, an 
appropriate threshold gray value for the preprocessed image 
is set. This results in a silhouette image. However, if the 
value is too high the profile becomes broken which prevent 
successful operation in the closing process discusseqmaemm 
Ghaptenr— 2. as shown in Figure 4.42. On the other hand 
decreasing the threshold value results in erroneous Prof ie 
as shown in Figure 4.43. In some cases, when the closing 
process is used, certain vital information is lost.  “Diggie 


effect can be seen in comparing the image in Figure *.44, 


LOC 





| 





and Figure 4.45. Therefore, loss of informations due to the 
attempt of eliminating noise and the inability in the 
preprocessing process to cope with the residual noise, yield 
erroneous profile. Consequently, maT ime woeecuns ine the 


succeeding classification steps. 
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Figure 4.42 Hargiee linings sine id. 





Figure 4.43 Low Threshold. 
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Figure 4.44 Before Closing. 





Figure 4.45 After Closing. 
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V. CONCLUSION 


The shape of the superstructure of a ship is the most 
important feature used in the classification algorithms. To 
extract the edge profile, a Sobel operator is employed. The 
limitation of a Sobel operator is that some small details of 
the edge is lost. For example, the image of a DD at a range 
of 79000 feet after applying the Sobel operator shows the 
superstructure and the radar. But the edge of a small mast 
disappeared as shown in Figure 5.1. Furthermore, if the 
threshold value is set too big, the edge image of the super- 
structure profile becomes broken. The top superstructure 
profile is obtained by setting the gray value under the 
slope between the bow and the stern to zero. We need to 
apply a contour tracking process to refine the superstruc- 
ture profile. For some images, the connection of the broken 
profile pieces may be achieved in a Closing operation as 
discussed in Chapter 2. The disadvantage of the Closing 
operation is that some small details of the profile may 
disappear. The superstructure profile of a freighter at a 
range of 53000 feet is shown in Figure 5.2. After the 
Closing process the detail of the cranes disappeared as 
Sfewm In Figure 5.3. 

The ship classification can be achieved by either the 
Fourier Coefficient method or the B-spline Coefficient 
method. In using these coefficients to classify ships, it 
is found that only the initial coefficients that lie between 
mee Oc to the 20th, are relevant, while the rest are not. 
Inspection of the comparison curves of the same class of 
ships shows that, similarities in patterns exists up to the 
2ZO-th point. Beyond that, diversities in shape are so great 


miac inclusion of those additional points for classification 
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will be of no use. Examples are shown in Figure 5.4 and 


Figuireworuo. 





Figures: 1 Edge Image of a DD at a Range of 79000 feet. 





Bacyu reo. 2 Before Closing Process. 
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Peegure i5.6 Atter Closime@e .6eess. 
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Logarithmic Magnitude of a CGN at a Range of 45000 feet. 


Figure 5.4 


LOG 


Logarithmic Magnitude of a CGN at a Range of 55000 feet. 


Figure 5.5 


In the B-spline coefficient method the technique of 
uneven knot selection has been employed. With the original 
Sampling points on the ship profile as input, a smaller set 
of approximation samples are collected. These can be used to 
reconstruct the ship profile with sufficient information 
retained from the original image. The execution time as 


approached to that of handling the original data directly 


has been greatly reduced. The reduction of the number of 
Sample points by almost a factor of 10 1s common. For a CGN 
Sameea set of original sampling points of 290, has been 


reduced to 36. 

Comparing the two methods, Fourler Coefficient and 
B-spline Coefficient, Ene sLormer method, in some cases 1s 
not effective to establish satisfactory classification of 
ships. This is due to difficulties in matching similarities 
of the shape of the coefficient curves. The latter, 
Meee) SULDasses the former in that it is able to classify 
Memes ships accurately using computer programs. Le is 
possible to improve the reliability of those two methods by 
Peducing noise in the data collection process and the 


preprocessing process. 
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APPENDIX A 
THE PROGRAM TO OBTAIN THE SUPERSTRUCTURE 
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6-Dec-19846 17:18:40 VAXK=- 
Source Listing 6-Dec-1984 17:218:31 _DRA 


program cutCinput,output»zainfile,outfil?a); 


type 
Byte = 0..255; 
imagerowl packed array [0..257] of bytes 
iumag2row2 packed array €0..255]) of anteager: 
1mag?row3 packed array (0..255] of byte; 
rowl = oacked serray €0..255] of brtes 
row2 packed array (0..255] of integer; 


var 
sopel : array (€0..63] of amagerow?2; 
is 3 : bytes 
f <: array C0..65]) of imagerowl1; 
outfile : file of amagerow3; 
Gusdysrangersdrightsnem sintegers 
infile : file of rowl; 
Slope : real: 
image ‘:array (€0..63] of rowl:; 
KeKlexlsyley2 = anteger; 
NUM1 ,NUM2 »NUM3»NUMS,maxesmingmaxl 3; ante2ger: 
cal sarray C€0..63] of rowl; 
thes : integer: 
NAME =: PACKED ARRAY [€1..20] OF CHAR; 


BEGIN 


HWRITELNC°INPUT SHIP FILENAME OUT CUT.DAT°)S 
READLNCNAMS)S 
WRITELNC “THRESHOLO”)>; 
READCTHES); 
WRITELNC°NUM TO CUT FROM LEFT”); 
REAQOCNUM1): 
WRITELNC°NUM2 TO CUT FROM RIGHT” )3 
REAOCNUM2Z); 
WRITELNC?NUM3 TO CUT FROM TOP”); 
READCNUM3)3 
WRITELNC°NUMS TO CUT FROM BUTTOM*): 
READ CNUMS); 
open CinfilesNAME, history ‘s=old, 
access_method ‘:=sequential, 
record_length :=256ysrecord_type s=fixed); 
open Coutfile, °CUT.dat”,mhastory :=newsrecord_lenaqth :=256, 
record_type :=fixed); 
resetCinfiled; 
resrite Coutfiled); 
as 0 
while not eof Cinfile) do 
oegin 
read Cinfile,image(Cij): 
for j3:=0 to 255 do 
fCi¢l,je¢l12 2= amageli, jJ3 
1:=i¢l1; 
end; 


Cfcompute sovoel*) 
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6-Ocec=196% 17275240 VAX- 
Source Listing 6~-Dec-1984% 17°18: 31 _ORE 


FOR J¢=0° 197 80 
FC64, Se 248 JS=FC635I+246835; 
for i s= 0 to 63 Go begin 
f€i*l,02 s= fli*l,1]9;5 
€Cit*el,257) 2:2££141,256)3 
ends 
for j:= 0 to 255 do begin 


{EO p37 14 >= f€1,j34139; 
f£C65,)41) s= F064, j+135 
end. 


109,03 -.=fC1,1 3, 

TCOs 257 JSS fh ls 256s: 
f065,032=f064,1]95 
£065 2257TIJ2=f£664,256)5 


C#set initial max,» min *) 


for i:= 0 to 63 do began 
for j:= 0 to 255 do began 

Oxz=fCi,» j3*+2=fCiel, j3*+f01+2,j33-fC1i,j+2) 
~m2efCielsj+2]~f0Ci+2¢ j*233 

dys=f{Cite¢2,jJ+2stCir+2,g+lJ+flit2, j+23-1 6145) 
~2e¢fCisej+lI-fCi»,j+2]9; 

sobelfisjJi=round( Cdxse2edyee2 )%20.5); 

if max < sobelfi,»zjj] then 


max := sobelfi,jiJ; 
if min > sobelfi»jJ then 
min := sobelli», jj; 
ends 
ends 


range = max-mins 
Cx rescalet) 


for 1°:= 0 to 63 do begin 
for j:= 0 to 255 do begin 
callfCi,sj3 *= round(CCsobelli,s, jJ—-min)*+255)/range), 
if €j3<=NUM1) or Cj>=255-NUM2) then 
calfi,j23 :=03 
if €i<¢=NUM3) of (€1>=63-NUMS) then 
calli»sj3 s= OF 


C2 cut theshold*) 
IF THES<>0 THEN BEGIN 
1f calfli,jjJ <=thes 


then 
callCisj} := 0 
else 
eallissd <= 255; 
END, 
ends 
ends 


Ct&find point at raw and aft*) 
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€=Hee-1984 7517-18250 VAX- 
Source Listing 6-Uec-137646 17:18:21 _ORA 


25 ors 

0 to 150 do begin 

or 1:=0 to 63 do begin 
r1f€calli,jJ=bright) and (n=0) then 


re © CO 
ee es 28 of 
“eo 


~os 3 
Ona 

—a- 3 . ce ee 
0 

re 2 


begin 
yl 25 je 
x1 c= is 
nm s= nels 


end; (tif =) 
if€calli,255—jj=bright) and (m=0) then 


begin 
Vico t= 2 39>) 5 
Mes lis 
ms= mel; 
ends; (8&1 f=) 


end; (=for) 
end; (Csfor=) 


BGerteln€~xi="exly”°yl= "syle, "x2=° ox2Zy9"yZ=°sy2)3 


C# cut lane) 


Slope :=1.0; 
af xl=x2 then 


begin 
slope <= 03 
for i:t=x1l+¢1 to 63 do began 
for js=yl-5 to y2+¢5 do 
calli,jJ :=0;3 
end:;(=for=) 


end; Czif3) 
if slope<>0 then 
begin 
Slope :=(x2-xl)/Cy2-yl): 
1f slope < 0 then begin 
slope := absCslope); 
for j := yl to y2 do begin 
x S= xl¢l-round(( j-yl)*slope); 
for i :=x to 63 do 
calli, jjJ:=05 
end.C=for=) 
end,.(+1f*) 
if slope >0 then begin 
for j:= yl to y2 do begin 
x 32 xnleleround(( j-~yl)*slope): 
for 1 :=x to 63 do 
calfa,yjJ 2:=03 
end:(+for=) 
end s( if) 
writeln(’step 2”); 
end. Ctifs) 


C#put outtiles) 


for 1:°:=0 to 63 do begin 


alae 


Source. Uis ting 


for 5 ==O0.%0 255 de 
outfiie*()] 32 calkes,3): 
put Coutfiled); 
end: (*forz) 
closaCianfile): 
end. 
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APPENDIX B 
THe PROGCKAM OF CONTOUR FOLLOWING 


13 





10=Bbec—-1 364 9)3..545-55 
Source Listang é-Dec-1986 17:38:13 


PROGRAM TRACECINPUT,OUTPUT sINFILE,OUTFILE)D: 


TYPE 
OYTEs= 205-255; 
IMAGEROW1 = PACKED ARRAY (£0..255]) OF BYTE: 
ROW = PACKED ARRAY (0..257) DF BYTE: 
VAR 
ReCeoFOeFleF2,F3eFSsF5eFb,FT,FBR => BYTES 
F s ARRAY (02-653 OF ROW: 
AsIMAGE 3 ARRAY £€0..633 OF IMAGERDWI; 
INFILE 3: FItLe OF IMAGERDW1I: 
R1,C1i1,R2,C2,cCOIF,I»,J,COUNT =: INTEGERS 
ROW,COL = ARRAY £0..512) OF INTEGERS 
DUTFILE =: FILE DE IMAGEROW1: 
NAME = PACKEO ARRAY [—1.--20) OF CHAR; 


PROCEDURE STORE: 


BEGIN 
CORETA:=C=1; 
ROACIJ>S=R-1; 
WRATELNC°COL=° ,COLTIIJ, “ROw=",ROWCII)D: 
COUNT <:=13 
Ts=I+1;3 

ENO; 


PROCEDURE CMOVE; 


BEGIN 
FQ SssFOCR—-1,C)3 FISHFCKR~-1.,C*1l13: F2SSFOR,»C+lIJs F3ISSFCR*+1,C4+1); 
FOS=FOR*1,CI3 FS2=FCR+1,C-1)3 F6S=FCR,C-133 FITZ=FCR-1,C-13;3 
Case COIR OF 
0: BEGIN 
IF €F3=0) AND CF2=255) THEN BEGIN C22Ce#13 CDOIR2S=135 END 
ELSE IF CF2=0) AND CF1=255) THEN BEGIN Ri=R-15 
STOREs Cs=zCels COIRS=03 ENO 
ELSE IF CF1=0) AND CFO=255) THEN BEGIN R:& 
COIR:=0%: END 
ELSE IF CFO=0) ANG CFT7=255) THEN BEGIN C:=C-15 
STORES Ri=R-15 COIR2=33 ONO 
ELSE IF CF7=0) ANDO CF6=255) THEN BEGIN C: 
COIRS=3: ENO 
ELSE IF CF6=C) AND CF5=255) THEN BEGIN Ri=Rel5 
STORES CssCH-1s COIR =33 END 
ELSE BEGIN Ri=Rel1l3 COIR2=23 END: 
ENO? 
ls BEGIN 
IF CF5=0) AND CF4=255) THEN 6GEGIN Ri=R15 
COIRS=2:3 END 
ELSE IF CF4=0) AND CF3=255) THEN BEGIN C2=Cel; 
STOREs R2=R#1i COIR =13 ENO 
ELSE IF CF3=0) ANC CF2=255) THEN BEGIN C2=Cel; 
COIR:=1: END 
ELSE IF CF2=0) ANC CF1=255) THEN BEGIN Ri=R—-13 
STORES s=C el.  -CGIRS=05) (END 
SLSE IF CFL=0) AND CF0=255) THEN BEGIN Ri=R-13 
COIR:=03. END 


if 


R-1¢ 


iT 
ao 

! 
" 
ee 


114 


VAX-11 
_ORAO: 


Source Lasting 


EUSE ITF CEQ=C) AND CF7=255) THEN 
STORES R2=R-1;s COIR2S=03 END 
ELSE BEGIN Cs=C-13 COIR? =33: ENO: 


ENO: 


23 


BEGIN 


BEGIN 


10-Dec-1984 13:53:33 


6-Lec-1968% 


Graco ts 


IF CFT7=0) AND CF6=255) THEN BEGIN C22C-1; 


COIR?S=3; ENC 
ELSE IF CF6=0) ANC CFS5=255) THEN 
STORE, C2=C-1¢5 COIRS=3; ENO 
ELSE IF CF5=0) ANO CFS=255) THEN 
COIR:=2%3 ENO 
ELSE IF CFS=0) ANC CF3=255) THEN 
STORE;s K2=Re1s COIR:=13 END 
ELSE IF CF3=0) ANG CF2=255) THEN 


BEGIN 


BEGIN 


BEGIN 


BEGIN 


REGIN 


BEGIN 


RS=Rel, 


Cls=c-is 


COIRs=13 END 
ESE 1F €F2=0) ANC CF1l=255) THEN 
STORE; C2=C+#1. COIR:=1;5 END 
ELSE BEGIN R2=R-15 COIRS=03s ENDS 
ENO; 
3: BEGIN 
IF CF1l=0) AND CFO=255) THEN BEGIN R2=R-1;% 
COIRS=0; ENO 
ELSE IF CFO=0) ANC CF7=255) THEN 
STORE, s=P—-15, COIR2=03; END 


ELSE IF CF7=0) ANC CF6=255) THEN 
COIR:=3¢; END 

ELSE IF CF6=0) ANG CF5=255) THEN 
STOREsS Ct=CHl1l;s COIR?=23 END 

ELSE IF CF5=0) ANC CF4=255) THEN 
EDIR2=23 ENDO 


ELSE IF CF4S&=0) ANC (CF3=255) THEN 
STORES Ri=Rel1, COIR? =2; END 
ELSE BEGIN Ci=C#15 COIRS=13; ENO; 
EN; 
END; 
ENDO; 


PROCEOURE INITIAL; 


BEGIN 
Por 1 S20 TO 255 O00 8SEGIN 
SGOLEIJ 3:=0; 
ROWCIJ °2=0;3 
ENDS 
FOR JjJ2=0 TO 63 OO BEGIN 
FOR 1:=0 10 255 OO 
ACR,CI:2=03 
ENDS 
es = 05 
END: 


PROCEOURE FIRSTs 


VAR 
No 
BEGIN 


M SINTEGER; 


BEGIN 
BEGIN 
BEGIN 


BEGIN 


LS) 


C2=c =); 
Rs=Rel135 
Rs=R+15 


C:=C+l; 


Li. 38513 


Vax%—11 
~OREOQ: 


L0-Cet-195S ae. 33 VAX~-11 


Source Lastang 6-Dec-198% 17238:13 _ORAD? 
N 2= OF M2=05 
FOR C :=Q0 TO 200 OO BEGIN 
FOR R :=0 YO 63 OO BEGIN 
IF CFER.~CIJ=255) ANO CN=0) THEN BSGIN 
Ci t= C$ Ri s=R5 NN S= Noi; 
ENO; 
IF CFOCR,255-C3=255) ANO (CM=0) THEN BEGIN 
C2 := 255-Cs R2 2= RE M 2EMOl1 SG 
ENO; 
ENOs 
ENO; 
WRITELNC °COL1=",Cl,°ROWL=°2R1, “COL2=°,C2, “ROW2Z=" »R2); 
ENO; 


CHSSVSE RAE AE SH SESS SESARKS ESE ARIAS ASS 
C# MAIN PROGRAM x) 
SEGIN 
WRITELNC°INPUT CUT OR CONLINE FILENANE”™) 3 
REAODLNCNAME)3 
OPEN CINFILEsNAME,HISTORY = OLO, 
ACCESS _ METHOO :=SEQUENTIAL, 
RECORO_LENGTH 2s=256,RECORO_TYPE S=FIXED): 
OPEN COUTFILE, “TR.~OAT” »HISTORY 2=NEW,RCCORO_LENGTH :=256, 
RECORGO TYPE <=SFIXEO); 
RESETCINGEILED; 
REWRITE COUTFILED: 
R :=0;% 
WHILE NOT EOF CINFILE) OO 
BEGIN 
READ CINFILE.-IMAGEECRI): 
FOR C := Q TO 255 DOO 
FOR#1,C413 2= IMAGECR,CI3 8 
R 3= Rl; 
ENO: 
FIRSTS 
KR S=R15. 6 3=€il. COIR 2-1. 
INITIAL? 
WHILE C€C<>C2) OO BEGIN 
STORES: 
CMOVEs; 
ENOS 
FOR Y:=0 TO COUNT OO 
ACRGWETI, COLE IM) 22255; 
FOR R:=0 T9 63 OO BEGIN 
FoR C:=C TO 255 QD 
OUTFILE*CCJ s= ACR,CI3 
PUTCOUTFILE)D: 
ENO? 
WRITELNC °NUMBER=°",COUNT) 3 
CLOSECINEFILE)DS 
ENO. 


6 


APPENDIX C 
THE SUBROUTINE "PARAM" 


analy 


NANDADNDADDNNDADANADNADAANNDNDADADANDAAANDANDADADDANADAAANAAAAAADAANAADANDADANAAANAAANAAH 


SUBRCUTINE PARAMC Xo VoeZowoMeZBeZE eK eSoNeoToCxXeCYeoFPeIOPT,IPAR, IER) 
GIVEN THE SET OF DATA POINTS CXCI)sYCII) WITH CORRESPONDING Z- 
VALUES ZCI) el=leZ2ee eM AND GIVEN ALSO THE SET CF PUSTIIVE 
NUMSERS wWCI) sl=le2veeeMe SUBROUTINE PARAM FINOS A SMOOTH APPROXIMAT-~- 
ING CURVE WITH PARAMETER REPRESENTATION X = SXCZ)- Y = SYCZ).~ 
S$x€Z> ANO SYCZ) ARE TWO SPLINE FUNCTIONS OF GEGREE K WITH THE NUMBER 
ANC THE POSITION OF THE KNOTS TC J) eJ=leo2eee-N AUTOMATICALLY 
CHOSEN SY THE ROUTINE. THE SMOOTHNESS OF SXCZ) AND SYCZ) IS 
ACHIEVED BY MINIMALIZING THE SUMCOXCRIRZS24¢0YCR) S42) WHERE OXCR) 
ANO OYCR) STAND FOR THE OISCONTINUITYJUMP OF THE KTH OERIVATIVE 
OF SXCZ) ANDO SYCZ) AT THE KNOT TCR), R=K42,---.N-K-1l. 
THE AMOUNT OF SMOOTHNESS IS OETERMINEO BY THE CONOITION THAT FCP) = 
SUMCWCTISCOXCID-SXCZC1))) 822 *#CYCII-SYCZCII)93%2)) BE <= Se WITH 
S A GIVEN NON-NEGATIVE CONSTANT. 
THE SPLINE FUNCTIONS SXCZ) AND SYCZ) ARE GIVEN IN THEIR B-SPLINE 
REPRESENTATION CB-SPLINE COEFFICIENTS CXCJ),RESP. CYC I) sJzleo we N-K-1) 
ANO CAN BE EVALUATEO BY MEANS OF FUNCTION OERIY. 
CALLING SEQUENCE: 

CALL PARAMCX gVeZoweMeZ5pZE aK oe SoNoleCXeCY,FPe,IGPT,IPAR,IER) 


INPUT PARAMETERS: 

x ARRAY,LENGTH My CONTAINING THE ASSCISSAE OF THE OATA POINTS 

ARRAY,LENGTH My CONTAINING THE ORDINATES OF THE DATA POINTS 
QARRAYMINIMUM LENGTH M,CONTAINING THE WEIGHTS WCI). 
INTEGER VALUE,CONTAINING THE NUMBER OF OATA POINTS. 
INTEGER VALUE,sCONTAINING THE OEGREE OF SXCZ) ANDO SYCZ).~ 
REAL VALUE,CONTAINING THE SMOOTHING FACTOR. 

INTEGER VALUE WHICH TAKES THE VALUE 0 OR 1. 

C: THE ROUTINE WILL RESTART ALL COMPUTATIONS. 

1: THE ROUTINE WILL START WITH THE KNOTS FOUNO AT THE 
LAST CALL OF THE ROUTINE. IF LOPT=1 THE DUTPUT 
PARAMETERS T ANO N ARE INPUT PARAMETERS AS WELL. 

IF JOPT=1 THE USER MUST PROVIOE WITH A COMMON BLOCK 
COMMON/OPTI/NROATACNEST) »FPO,FPOLOeNPLUS 
IPAR = INTEGER FLAG. 


Y 
id 
M 
K 
S 
IOPT 


IOPT= 
TORT= 


IPAR = QO: FOR EACH DATA POINT CXCT)+YCID) THE PROGRAM AUTIMATICALLY 


CHOCSES A CORRESPONOING VALUE OF THE PARAMETER Zy I.E 


ZCLIFOSZCII=AZCIA“1LI+SQRTICCXCI “XC I-19) Sa2e CVC II-V¥C1-1)) 452) 
THE BOUNDARIES FOR THE PARAMETER Z ARE CHOSEN AS FOLLOWS 


78 = 201) 5 Z2ZE = 2040. 
IPAR = 1: THE USER HIMSELF PROVIOES WITH THE VALUES OF THE 
PARAMETCR Z ANO WITH THE BOUNOARTIES ZB ANO ZE. 
M > ARRAY,LENGTH Mey CONTAINING THE VALUES OF THE PARAMETER 2 
CIPAR = 1) 
Z8,ZE: REAL VALUES, CONTAINING THE BOUNOARIES OF THE PARAMETER 2 
CIPAR = 1). 


OUTPUT PARAMETERS: 


T > ARRAY,LENGTH NEST CSEE OATA INITIALIZATION STATEMENT), 
WHICH CONTAINS THE POSITION OF THE KNOTS»I.~E~ THE POSITION 
OF THE INTERIOR KNOTS TCK4+2),ceeTON-K-1), AS WELL AS THE 
POSITION OF THE KNOTS TC1)=TC2)=.2.-=7CK*1)=72B8 AND ZE = 
TON-K)=.ee=TCN) WHICH ARE NEEOEO FOR THE B-SPLINE REPRESENT. 
CXe,CY: ARRAYS,»LENGTH NEST, CONTAINING THE B-SPLINE COEFFICIENTSS 
OF SXCZ)- RESP. SYCZ).~ 
INTEGER VALUE,CONTAINING THE TOTAL NUMSER OF KNOTS. 
REAL VALUE, WHICH CONTAINS THE SUMCWISCXI-SXCZI))%42) 
* SUMCWISCYIT-SYCZI))342) 5 L=1,29---M~ 
IER =: ERROR COOE 
IER=0: NORMAL RETURN. 
ITER=-12°NORMAL RETURN,SXCZ) AND SYCZ) ARE INTERPILATING SPLINES 


N 
FP 


@e &¢6@ 
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0930 
0940 
0950 
0960 
9970 
0980 
0990 
1000 
1010 
rez 
1030 
1040 
1050 
1060 
1070 
10380 
1090 
1100 
1110 
1120 
1130 
1140 
1150 
1160 
1170 
1180 
ji re 6) 
1200 
1210 
T2200 
yg 110 
1240 
iZou 
1260 
1270 
1280 
1290 
13900 
1310 
1320 
1330 
1340 
1350 
1360 
L370 
1380 
1390 
1400 
1410 
1420 
1430 
1440 
1450 
1460 
1470 
1480 
1490 
1500 
1510 
1520 
1330 


——— 


ANAAMAAANAARADADAANANAANDANDNAANAN 


NAAANDANNAANAAN 


AQOAAANDANANAANNAAN 


aan 


PER=—Z>NORMAL RETURN, SXCZ) ANO SYCZ) ARE POLYNOMIALS OF OEGREE K 


ITER>O ABNORMAL TERMINATION. 

DPER=LSTHE REQUIRED STORAGE SPACE EXCEEOS THE AVAILABLc 
STORAGE SPACE,SPECIFIED BY THE PARAMETER NEST. 
PROBABLY CAUSES:NEST OR S TOO SMALL. 

TER=254 THEORETICALLY IMPOSSIBLE RESULT WAS FOUND DURING 
THE ITERATION PROCESS. 

PROBABLY CAUSES: TOL TOO SMALL 

TER=32THE MAXIMAL NUMBER OF ITERATIONS MAXIT HAS BEEN 
REACHED. 

PROBABLY CAUSES:MAXIT OR TOL TOO SMALL.~ 

IER=10°SOME OF THE INPUT DATA ARE INVALIOCSEE RESTRICTIONS).~ 


RES URICTIONSS 
1) M> K > QO 
meee e <= LCR) < ZOER41) <= ZE, ReledeeweM—l1l. CIPAR = 1) 
3) WOR) D> Ov RE1,2,---H~ 
4) S >= 0. 
5) NEST D= 24K2. 


OTHER SUBROUTINES REQUIRED: 

BSPLIN,COSSIN»,ROTATE,BACKeNKNOT,OISCO ANO RATION. 

OIMENSION XCM) YOCOM) oeWOM) eo 2CM),TC 200) CXC 200)-CYC200), 
< FPINTC200)eRXC200) »RV¥C 200) 20ITAG C200) sOPRIMECZ00),» 

€ 66200 46).8C20007).QCS00 sb) eMC TI »>NROATAC2Z00),AC200.5) 

COMMOCN/ZOPTI/SNROATACNEST) sp FPO eFPOLD eNPLUS 

NROATA? INTEGER ARRAY,LENGTH NEST, WHICH GIVES THE NUMBER OF 
OATA POINTS INSIO€G EACH KNOT INTERVAL. 

FPQ >: REAL VALUE, WHICH CONTAINS THE SUMCWISCXI-SXCZ1I)) S42)+ 
SUMCWISCYI~SYCZID)I%22) WITH SX€Z) ANO SYCZ) LEAST-SQUARES 
POLYNOMIALS OF OEGREE K. 

FPOLO =: REAL VALUE.WHICH CONTAINS THE SUMCWIECXI-SXC2Z1))%52)+ 
SUMCWISCYI-~SYC2Zi))4%2) WITH SXCZ) AND SYCZ) LEAST-SQUARES 
SPLINE FUNCTIONS CORRESPONOING TO THE LAST FOUNO SET OF 
KNOTS Sut ONE. 

NPLUS =: INTEGER VALUE,GIVING THE NUMBER OF KNOTS OF THE LAST 
SET MINUS THE NUMBER OF THE LAST SET BUT ONE. 

COMMON/OPTI/SNROATA,FPO,FPOLD,NPLUS 
OATA INITIALIZATION STATEMENT TO SPECIFY 

TOL : THe REQUESTED RELATIVE ACCURACY FOR THE ROOT OF FCP) = §. 

MAXIT:S THE MAXIMAL NUMBER OF ITERATIONS ALLOWEO. 

NEST = AN OVER-ESTIMATE OF THE NUMBER OF KNOTS Nj THIS PARAMETER 
MUST BE SET B8Y THE USER TO INOICATE THE STORAGE SPACE 
AVAILABLE TO THE SUBROUTINE. THE OIMENSION SPECIFICATIONS 
OF THE ARRAYS TeCXeCYeNROATA,FPINT,RXRY,OIAGsOPRIMECN) ,g 
ACN eK) sGONe Kl) eBCNeK*2) -CCMeX*1) ANDO HCK*2) OEPENO 
ON NeM ANO Ke SINCE N IS UNKNOWN AT THE TIME THE 
USER SETS UP THE OIMENSION INFORMATION AN OVER-ESTIMATE 
OF THESE ARRAYS WILL GENERALLY BE MADE. THE FOLLOWING 
REMARKS ARE INTENOEO TO HELP THE USER 

1) 26K+2 <= N C= MeKel] 
2) THE SMALLER THE VALUE OF S, THE GREATER N WILL BE. 
3) NORMALLY N = M/2 IS AN OVER-ESTIMATE. 
OATA TOLJO.OOL/ »pMAXITS20/ gNEST/S200/ 
BEFORE STARTING COMPUTATIONS A DATA CHECK I35 MADE. IF THE INPUT 


DATA ARE INVALIO CONTROLE IS IMMEDIATELY REPASSEO TO THE ORIVER 
PROGRAM CIER=10). 


IER = 0 
Kl = Ke] 
K2 = Kile] 


NMIN = 22K] 


LES 


1540 
1550 
1560 
i570 
1580 
L590 
1600 
1610 
1620 
1630 
1640 
1650 
1660 
1670 
1680 
1690 
1700 
1710 
L720 
PT 30 
1740 
Li oe 


1790 
1800 
1810 
1820 
1830 
1840 
1650 
T1860 
1870 
reso 
1830 
1900 
1910 
LaZ0 
1930 
1940 
1950 
1960 
1970 
1980 
HOSTS 
2000 
2010 
Z0ce 
2030 
2040 
205.0 
2060 


2080 
2090 
2100 
2110 
2120 
2130 
214C 


TECKR.CES02 LER ae 
IFCM.LT-K1 .OR. NEST.LY<NAMIN)D JER = 20 
TFCS.k712.Ge) TER = FF 
IFCIER.NE.~0) GO TO 440 
C CHECK WHETHER THE Z-VALUES ARE PROVIOEO WITH BY THE USER. 
IFCIPAR.NE.O) GO TD 6 
C FINO FOR EACH OATA POINT A CORRESPONOING VALUE OF THE PARAMETER Z 
C AND FIX THE BOUNDARIES ZB ANO ZE.j 
ZG1D> = 10. 
oo & I=2, 
ZCLT) = ZCIKLI*SQRICCXCIINXCI~1) S424 CVC IT -YCI~-1))*32) 
4 CONTINUE 
ZB = 71) 
ZE = 2C() 
6 IFCZE.GT.ZC1) ~CRe ZESLIAZCM) sORS WO DSL. 0) eee 
00 10 I=2,M 
IFCZCI=-1L).GE<ZC1) .ORa. WCIDSEESCOS) TER = 920 
10 CONTINUE 
IFCIER.NE.0) GO TO 440 
C CALCULATION OF ACC, THE ABSOLUTE TOLERANCE FOR THE ROIJT OF FCP)=S. 
ACG = TOESS : 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCECECEC EG CC ECCCECLCCEC CC CCL ee eC eee CeCe Creme aes 
PART 1: OETERMINATION OF THE NUMBER OF KNOTS ANDO THEIR POSITION C 
SVEVBERPSSSE AK AESSEZESSASESA SK AIESSEVA SASHA SESHSESAHEK RH AKHESSESESTHAREER C 
GIVEN A SET OF KNOTS wE COMPUTE THE LEAST~SQUARES SPLINES SXINFCZ) C 
ANO SYINFCZ).IF THE SUM FCPSINF)<C=S we ACCEPT THE CHOICE OF KNOTS. c 
OTHERWISE WE HAVE TO INCREASE THEIR NUMBER. C 
THE INITIAL CHOICE OF KNOTS OEPENDS ON THE VALUE OF S AND IOPT. G 
IF S=0 WE HAVE SPLINE INTERPOLATION: IN THAT CASE THE NUMBER OF G 
KNOTS EQUALS NMAX = MeKely G 
IF § > 0 AND C 
ITOPT=0O WE FIRST COMPUTE THE LEAST~SQUARES POLYNOMIALS OF C 
CEGREE K3 N = NMIN = 2€K2 C 
TOPT=1 WE START WITH THE SET OF KNOTS FOUNO AT THE LAST C 
CALL OF THE ROUTINE, EXCEPT FOR THE CASE THAT S > FPO: THEN c 
WE COMPUTES OIRECTLY THE LEAST-SQUARES POLYNOMIALS OF OEGREE KK. C 
CCCCCCCCCCCCCCCCCCCCCCCCECCCCECCCECECE COCEC CCE CCE CeCe eee ere errr Greece eee 
DETERMINE NMAX, THE NUMSER OF KNOTS FOR SPLINE INTERPOLATION. 
NMAX = M4eKi 
TECS.«61.0.) GO 10745 
C IF S=0, SxXCZ) AND SYCZ) ARE INTERPOLATING SPLINES. 
N = NMAX 
C TEST WHETHER THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE ONE. 
IFCN.GI.-NEST) GO TO 420 
C FINO THE POSITION OF THE INTERIOR KNOTS IN CASE OF INTERPOLATION. 
MK1l = M~K] 
IFCMKI1.EQ.0) GO TO 60 


DOAANAAANAAANAANAAANANM 


K3 = K/2 
I = K2 
Js KS42 


IFCK3¥22.EC-K) GO TU 30 
O00 20 L=1,MK1 
TcI)D = ZCJ) 
] = i448 
J = 52} 
20 CONTINUE 
GO TO 60 
30 DO 460 L=1,MKi 
TCL) = C2 Gs - 20-1) ) 302 5 
I = Il 
J 


20 


2150 
2160 


21380 
2190 
2200 
2218 
as 74,0) 
2230 
2240 
2250 
2260 
2270 
2280 
2290 
2390 
2310 
2320 
2338 
2340 
2350 
2360 
2310 
2380 
2378 
2400 
2410 
2420 
2430 
2440 
2450 
2460 
2470 
2480 
2490 
2500 
2510 
2520 
2530 
2540 
2530 
2560 
2570 
2580 
253990 
2600 
2610 
2620 
2630 
2640 
2650 
2660 
2670 
2680 
2690 
2700 
zie 
2720 
2730 


2740 
2760 


60 CONTINUE 2160 


GO TO 60 270 

Cc IF S30 OUR INITIAL CHOICE OF KNOTS CEPENOS ON THE VALUES OF IOPT. 2790 
C IF TOPTH0 GR IOPT=A1 ANDO S>=FPO, WE START COMPUTING THE LEAST-SQUARES 2470 
c POLYNOMIALS OF OEGREE K WHICH ARE SPLINES WITHOUT INTERIOR KNOTS. 2800 
C IF TOPT=1 ANC FPODS WE START COMPUTING THE LEAST-SQUARES SPLINES 2810 
C ACCOROING TO THE SET OF KNOTS FOUNO AT THE LASY CALL OF THE ROUTINE. rays) 
Som 1FCIOPT.LE.0) GO TO $0 2830 
IFCFPO.GT.S) GO TO 60 2840 

SO N = NMIN Z2e50 
NROATAC1) = M~2 2860 

C MAIN LOOP FOR THE OIFFERENT SETS OF KNOTS.» MIS A SAVE UPPER BOUND 2870 
C FOR THE NUMBER OF TRIALS. 2880 
60 00 200 ITER = 1,5 2890 
IFCN-EQeNMIN) IER = -2 2900 

C FINO NRINT. TNE NUM3SER OF KNOT INTERVALS. 2910 
NRINT = N-NMINe@] 2920 

C FINO THE POSITION OF THE ADDITIONAL KNOTS WHICH ARE NEEDED FOR 2330 
C mae B8=-SPLINE REPRESENTATION DF SXCZ) ANDO SYCZ). 23940 
NKi = N-~-Kl 2950 

I = N 2960 

00 70 J=1,K%1 2970 

TOJD = 28 29 8:0 

held) = ZE 2930 

I = J-1 3000 

70 CONTINUE 3010 

c COMPUTE THE B-SPLINE COEFFICIENTS OF THe LEAST—SQUARES SPLINES SXINFCZ) 3020 
C ANG SYINFCZ)~. THE OBSERVATION MATRIX A IS BUILT UP ROW BY ROw ANDO 3030 
cE REOUCED TO UPPER TRIANGULAR FORM 3Y GIVENS TRANSFORMATIONS 3040 
C WITHOUT SQUARE ROOTS.j AT THE SAME TIME FP=FCP=INF) IS COMPUTED 3050 
aea= OO. 3060 

C INITIALIZE THE OBSERVATION MATRIX A. 3070 
09 890 I=1,NK1 3030 

DIAGCI)D = 0. 3090 

RXCI) = O. 3100 

RYCI) = 0. 3110 

0D 80 J=1leK 3120 

ACI,J) = 0. 3130 

80 CONTINUE 3140 

L = Al 3150 

00 130 IT=1,M 3160 

Ce FETCH THE CURRENT OATA POINT XCIT)eYCIT) eZCIT). 3170 
Lins XC PT) 3180 

YI = YCIT) 3190 

ZI = ZEIT) 3200 

WI = wCIT) 3210 

C SEARCH FOR KNOT INTERVAL TCL) <= ZI <= TCL41). 220 
iat leaGeeNCL#l]) sANO. LeNE.NKI1) L = Le] 3230 

c EVALUATE THE (€K#1) NON-ZERO B-SPLINES AT ZI ANO STORE THEM IN Q. 3240 
CALL BSPLINCTeNeKeZI pi ot) 3250 

00 90 [=1,K1 3260 
PECHCID~.LT~O.1E-07) HCI) = O. 3270 

QcIT,I) = HCI) 3280 

90 CONTINUE 3290 

c ROTATE THE NEw ROW OF THE OBSERVATION MATRIX INTC TRIANGLE BY 3300 
| C GIVENS TRANSFORMATIONS WITHOUT SQUARE ROOTS. 33710 
J = L~fl1 3320 

| DO 110 IT=1,€1 3330 
| IFCWI2EQ.0.2.) GO TO 130 3340 
| J = J+1 3350 
| PIV = HCI) 2340 


ah 





LFCPIV. 6020-0 GOs 1Gr110 3370 


C CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION. 3360 
CALL COSSINCPIV.WI,OIAGCI) »COS,SIN)D 3390 

C TRANSFORMATIONS TO RIGHT HANO SIOES. 3400 
CALL ROTATECPIV,COSsSIN, XI eRXCI)) 3410 

CALL RITATECPIV, COS «SIN, Yleek VCO 3420 

IFPCICEQ. KID GO 10-120 3430 

I2 = 0 3440 

I3 = Il 3450 

00 100 I1 = I3,K1 3460 

125] l2e 3470 

C TRANSFORMATIONS TO LEFT HANO SIDE. 3480 
CALL ROTATECPIV.COSsSINSHCI1) ACI ,12)) 3490 

100 CONTINUE 3500 
110 CONTINUE . 3510 
C AOO CONTRIBUTION OF THIS ROW TO THE SUM OF SQUARES OF RESIOUAL 3520 
C RIGHT HANO SIDES.~ 3530 
120 FP = FPeWIl aC xl s=2¢YIs22) 3540 
130 CONTINUE 3550 
TPCLEx.EGa-2) FROG= sre 3560 

C BACKWARO SUSSTITUTION TO OBTAIN THE B-SPLINE COEFFICIENTS. 3570 
CALL SACKCA,RX GNK Ihe C XD 3580 

CALL BACKCARY,NK1,K,CY) 3590 

C TEST WHETHER THE APPROXIMATION X=SXINFCZ) + Y=SYINFCZ) IS AN 3600 
C ACCEPTABLE SCLUTION. 3610 
FPMS = FP=-§ 3620 
IFCABSCFPMS).LT-.ACC) GO TO 440 3630 

C IF FCP=INF) < S ACCEPT THE CHOICE OF KNOTS. 3640 
LECFPMS. bl. 0. 8 GO bo 250 3650 

C IF N=NMAX,.SXINFCZ) AND SYINFCZ) ARE INTERPCLATING SPLINES. 3660 
IFCNLEQLNMAX) GO TI 430 3670 

C INCREASE THE NUMBER OF KNOTS. 3680 
C IF N=NEST WE CANNOT INCREASE THE NUMBER OF KNOTS BECAUSE OF 3690 
C THE STORAGE CAPACITY LIMITATION. 3700 
IFCNSEQ.NEST) GO TO 420 3710 

C DETERMINE THE NUMBER OF KNOTS NPLUS WE ARE GOING TO ADO. 3720 
IFCIER.EG.0) GO TO 140 3730 

NPLUS = 1 3740 

IER = 0 3750 

GO To 150 3760 

140 NPL1 = NPLUS#2 37710 
IFCFPOLD-FP.GT.~ACC) NPL1 = FLOATCNPLUS)SFPMS/CFPOLO-FP ) 3780 

NPLUS = MINOCNPLUS*2 »MAXOCNPL1,NPLUS/2,41)) 3790 

150 FPOLD = FP 3800 
C COMPUTE THE SUMCWISCCKI-SXINFCZIDS42*¢CVYI-SYINFCZI))2%*2)) FOR 3810 
C EACH KNOT INTERVAL TCJ*K) C= ZI C= TCS*#K41) ANO STORE IT IN 3820 
C ‘FRINTC J) eJHale2ee eo NRINIT.« 3830 
FPART = Q. 3840 

I = 1 3850 

L = K2 3860 

NEW = O 3B70 

00 180 IT=1.™ 3880 

IF CPCI «LI. TCU) «OR. LoGIONK? Gor 1Oece 3890 

NEW = 1 3900 

L = tl 3910 

160 TERM] = 0. 3920 
TERM2 = O. 3930 

LO = L-K2 3940 

00 170 J=1,.Kl 395C 

LO = LO) 3960 

TERM1 = TERMIL*CXCLO)SQlIT. J) TG 


NA 


C 


C 


C 


C 
C 


enrieee= VER MIeCY CEC osC Cl T. J) 


170 GENT Tive 


TERM = WCITISCCTERMI-XCIT))#S2+¢CTERM2-YCIT)) *%2) 
FPART = FPART+4TERM 

IFCNEWeEQ.0) GO TO 180 

STORE = TERM*0.5 

FPINTCI) = FRPART=STORE 

I = JI41 

FPART = STORE 

NEw = 0 


180 CONTINUE 


FPINTCNRINT) = FPART 
OO 190 L=1,NPLUS 


AODG & NEW KNOT. 
CALL NKNOTCZLeMeTeNeFPINT ,WROATA,NRINT ) 
TEST WHETHER WE CANNOT FURTHER INCREASE THE NUMSER OF KNOTS. 
IFCNLEQLNMAX LOR. NeEQ.NEST) GO TO 200 
190 CONTINUE 
RESTART THE COMPUTATIONS WITH THE NEW SET OF KNOTS. 


200 CONTINUE 
TEST WHETHER THE APPROXIMATION X=SKC2Z)eV=SYCZ) WITH SxCZ2) ANDO SYC2Z) 


THE LEAST=-SQUARES KTH OEGREE POLYNOMIALS. IS A SCLUTION OF OUR PROBLEM 
250 IFCIER.EQ--2) GO TO 440 


fummmmeminecccGcCCCCCCCCCCCECCCCCCCCCCCCCECCCCCECCECCCECCCCCCCCCCCCCCCCCCC 


AMNATNAANAANNAANANANAANANAAAANN 


- 


PART 2: OETERMINATION OF THE SMOOTHING SPLINES SXPC2) ANO SYPC2). 
SRAER KHER GK SESERK SSAA SRSA SSSA SE SESH SKS ASRAS AE SHAG SHAS AS SERS ALTE SS 
WE HAVE DETERMINEO THE NUMBER OF KNOTS AND THEIR PCSITION. 

Meeewicow COMPUTcE THE B-SPLINE COEFFICIENTS OF THE SMOOTHING SPLINES 


SxXxPCZ) AND SYPCZ). THE OBSERVATION MATRIX AI3 EXTENDED BY THE ROWS 


Oe MATRIX 8 EXPRESSING THAT THE KTH GERIVATIVE OISCONTINUITIES OF 
SXPCZ) AND SYPCZ) AT THE INTERIOR KNOTS TCK4*#2),--.-TON-K-1) MUST BE 
ZERO~} THE CORRESPONOING WEIGHTS OF THESE ADDITIONAL ROWS ARE SET 
TO I/SQORTCP). ITERATIVELY WE THEN HAVE TO DETERMINE THE VALUE OF P 
SUCH THAT FCP)=SUMCWI2CCXI-SXPCZIDDS424¢CVYI-SYPCZI))IFF2) BE = Se WE 
ALREAQY KNOW THAT THE LEAST-SQUARES POLYNOMIALS CORRESPOND TO P=0, 
AND THAT THe LEAST-SQUARES SPLINES CORRESPONO TO P=INFINITY. THE 
ITERATION PROCESS WHICH IS PROPOSED HERE, MAKES USE OF RATIONAL 
INTERPOLAYION. SINCE FCP) 15 A CONVEX ANDO STRICTLY DECREASING 
FUNCTION OF P, IT CAN BE APPROXIMATEO BY A RATIONAL FUNCTION 

ROP) = CUZP4V)/CP4W). THREE VALUES OF PCP1,.P2,.P3) WITH CORRESPOND- 
Mise VGCUES OF FCP) CFISFCP1)—-SeF22FCP2)-S,F3=FC93)-5) ARE USED 

TO CALCULATE THE NEw VALUE GF P SUCH THAT RCP)=S. CONVERGENCE IS 
GUARANTEEOQ BY TAKING F1>0 AND F3<O. 


M~pmmeeeocetececcCCCCCCCCCCCCCCCCCCCCCCCCCECCCCECCCCCCCCCCCCCCCCCCCCCCCCC 


EVALUATE THE OISCONTINUITY JUMP OF THE KTH DERIVATIVE OF THE 

B—-SPLINES AT THE KNOTS TCL) »L=K02,2...N—-K-1 AND STORE IN Bo 
CALL CISCOCT.N»K2,8) 

INITIAL VALUE FOR P. 


Brees Os 

FL = FPO-S 
P3 = -1. 
F3 = FPMS 
P= -FI/F3 
ICHECK = 0 


N8 = N-NMIN 
ITERATION PROCESS TO FINO THE ROOT OF FCP) = S. 
OO 350 ITER=1,MAXIT 
THE ROWS OF MATRIX 2 WITH WEISHT 1/S5SQRTCP) ARE ROTATED INTO 
THE TRIANGULARISEO OBSERVATION MATRIX A WHICH IS STORES IN Ge 
PINV = 1.-0/P 
OO 260 [=1.NK1 


1Z3 


\@ 
C 
c 
C 
C 
c 
C 
Cc 
6 
C 
C 
C 
C 
C 
C 
c 
: 
C 
‘G 
C 


3530 
39390 
4000 
4010 
4020 
4030 
4040 
4050 
4050 
4070 
40380 
4090 
4100 
4110 
4120 
4130 
4140 
4150 
4150 
4170 
4180 
4190 
« 200 
4210 
4220 
4230 
4240 
4250 
4260 
4270 
4280 
4290 
4300 
4310 
4320 
4330 
4340 
4350 
4360 
4370 
4380 
4390 
4490 
4410 
4420 
4430 
4440 
4450 
4460 
4470 
4480 
4490 
4500 
4510 
4520 
4530 
4540 
4550 
4560 
4570 
4ear 


UPRIMECI)D = OTAGCI) 4590 


Exc] ) = 72x) 4600 

CYCI) = RYCIS “610 

GCI sk1> = 0c 4620 

00 240 J=1.K 4630 

GcI,J) = ACI,J) 4640 

260 CONTINUE 4650 
96 300 11=1,N6 “660 

C THE ROW OF MATRIX B IS ROTATED INTO TRIANGLE BY GIVENS TRANSFORMATIONS. 4670 
GO 270) 1=1,2 “620 

HCl). = BCiig i) £690 

270 CONTINUE 4700 
x] =20. 4710 

Y] <=""0.. “720 

WI = PINV 4730 

00° 290 J=1 tend 4740 
IFCWILEQ.0.) GO TO 300 “750 

PIv = H(1) 4760 

C CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION. 4770 
CALL COSSINCPIVewWI,OPRIMECJ),COS,SIN) 4780 

C TRANSFORMATIONS TO RIGHT HANO SIDES. , 4790 
CALL ROTATECPIV,COSsSINAXICXC)) 4800 

CALL ROTATECPIV,COS,SIN,YI~CYCJ) ) 4810 

IFC J-EQ-NK1) GO TO 300 4820 

l20=9 4830 

IFCJ-GTNB) I2 = NK1I-J 4840 

00 280 I=1,12 4850 

C TRANSFORMATICNS TO LEFT HANO SIDE. 4840 
CALL ROTATECPIV,COS,SINSHCI+1) -GCJ,1)) 4870 

HCI) = HCI41) 4830 

280 CONTINUE 4890 
HCI241) = 0. “900 

290 CONTINUE 4910 
300 CONTINUE 4920 
C BACKWARO SUBSTITUTION TO OBTAIN THE 8-SPLINE COEFFICIENTS. 4930 
CALL BACK CG,CXeNK1,K% CX) 4940 

CALL BACKCGsCYoNK15K&,CY) 4950 

C COMPUTATION OF FCP). 4960 
Fe = Q. 4970 

L = K2 4980 

00 330 IT=1,.m 4990 
IFCZCIT)SLTATCL) .OR. LGTONK1) Gout 0310 5090 

tL = Ll 5010 

310 LO = L=-n2 5020 
TERM1 = 0. 5030 

TERM2 = 0. 5040 

00 320 J=1,K1 5050 

LO = LO) 5060 

TERML = TERMI4+CXCLOISQCIT, I) 5070 

TERM2 = TERM24CYC(LO)#QCIT,Y) 5030 

320 CONTINUE 5090 
FP = FPeW CIT SC CTERMI-XCIT))2824¢C(CTERM2-YCIT) #22) 5100 

330 CONTINUE 3110 
C TEST WHETHER THE APPROXIMATION X=SXPC(Z).Y=SYPC(Z) IS AN ACCEPTABLE 5120 
C SOLUTION. 5130 
FPMS = FP-S 3140 
IFCABSCFPMS).LT-ACC) GO TO 440 5150 

C TEST WHETHER THE MAXIMAL NUMBER OF ITERATIONS IS REACHED. 5160 
LPCITERsEO. MAKIN) GO to. 6o0m 5170 

C CARRY OUT ONE MORE STEP OF THE ITERATION PROCESS. £12380 
P2 = P 51990 


124 


F2 = FPMS 
PEGlCCHECK ONE. 0) GO -Td. 340 
TECCE2=-F3).GT.~aACG) GO TO 335 

Meeourm INITIAL CHOICE OF P IS TOO L&arGé. 
P = P#0Q.1E-02 
P3 = P2 
=3 = F2 
GO TO 350 

335 IFCCF1I-F2)-GT.aCC) GO YO 340 

Cy OUR INITIAL CHOICE OF P XS TOO SMALL 
TYPE *,°VALUE OF P°,P 
P = P*J.1E*04 


Pl = P2 
Fl = F2 
GO TO 350 


C TEST WHETHER THE ITERATION PROCESS PROCEEOS AS THEORETICALLY 
C EXPECTED. 
340 TRC? cGESFL ~OR. FZ2.LE.F3) GO TO 4610 
ICHECK = 1 
C FINO THE NEW VALUE FOR P. 
P = RATIONCP1,F1l,P2,F2,P3,F3) 
350 CONTINUE 
c ERROR COOES ANO MESSAGES. 
400 IER = 3 


GO TO 4640 
410 IER = 2 
GO TO 440 
420 IER = 1 
GO TO 440 
430 JER = -] 
440 RETURN 
ENO 


SUSROUTINE BSPLINCT se NoKe Xot oH) 
SUBROUTINE BSPLIN EVALUATES THE (CkK#1) NON-ZERO B-SPLINES OF 
DEGREE K AT TCL) C= X K€ TEL#1) USING THE STABLE RECURRENCE 
RELATION OF OE S8OOR AND COX. 
THE OIMENSTION SPECIFICATIONS OF THE FOLLOWING ARRAYS MUST BE 
AT LEAST HC Kel) eHHCK). 
OIMENSION TCN).HC6)eHHCS) 
aC1) = 1. 
00 20 J=1,X* 
00 10 IT=1l,J 
HHCI) = HCI) 
10 CONTINUE 
HCl) = 6. 
QGe20 l=1.5 J 
LI = tel 
LJ = LI-J 
F = HHCID/CTICLID-TCLI)) 
HCI) = HCI)+F eC TCLID-X) 
MCI*¢1) = FSCx-TCLJ)) 
20 CONTINUE 
RETURN 
ENO 
SUBROUTINE COSSINCPIV,WI»WwWsCOS.SIN)D 
c SUBROUTINE COSSIN CALCULATES THE PARAMETERS OF A GIVENS 
C TRANSFORMATION WITHOUT SQUARE ROOTS. 
STORE = PIVawWI 
00 = Ww+STORE*PIV 
IFC ABSCOOD.LT-1-£&-36) OO=1.E-36 
CCS = wws0d 


NAANN 


ZS 


On AaAaANY 


oan 


SIN = STOREZOC 


WwW = DOD 

WI = COS=wI 
RETURN 

END 


SUBROUTINE ROTATECPIV.,COS,SIN,»A,3S) 


SUBROUTINE ROTATE APPLIES A GIVENS ROTATION TO A ANC 8B. 


STORE = 8 

B = COS*STORE*SINEA 

A = A-PIVESTORE 

RETURN 

ENO 

SUSROUTINE BACKCA,2 NX 5C) 


SUBROUTINE SACK CALCULATES THE SOLUTION OF THE SYSTEM OF 
EQUATIONS A=C = Z WITH A AN X N UNIT UPPER TRIANGULAR MATRIX 


OF SBANOWIOTH Kel. 


ATTENTION? THE FIRST DIMENSION SPECIFICATION OF MATRIX & MUST 


BE THE SAME &AS IN THE CALLING PROGRAM, 
DIMENSION AC200,K),Z20N) CCN) 
CCN) = ZCN) 
I = N-1 
IFCI.EQ.0) GD TO 30 
DO 20 J=25N 
STORE = ZCI) 
Il = K 
LEC I.LE CKD elise 2-1 
M= I 
00° 10 Wate) 
M= M+} 
STORE = STORE-CCM)ISACI gL) 
10 CONTINUE 
CCI) = STORE 
I = fei 
20 CONTINUE 
30 RETURN 
END 
SUSRDUTINE NKNOTCX Me TeNeoFPINT»sNROATA,NEINT ) 
SUBROUTINE NKNOT LOCATES AN ADDITIONAL KNOT FOR 
K AND ADJUSTS THE CORRESPONDING PARAMETERS,..£.~ 


A SPLINE DF OcGREE 


T - CHE, POSITION OF THE KNOTS. 

N > THE NUMBER OF KNOTS. 

NRINT 3 THE NUMBER OF KNOTINTERVALS. 

FPINT > THE SUM ODF SQUARES OF RESIOQUAL RIGHT HAND SIDES 


FOR EACH KNOT INTERVAL. 


NROATAS THE NUMBER OF DATA POINTS INSIOE EACH KNOT INTERVAL. 
THE ARRAYS T,FPINT ANO NROATA MUST HAVE THE SAME DIMENSION 


SPECIFICATIONS AS IN THE CALLING PPOGRAM. 
DIMENSION X€M),TC200),FPINTC200) sNROATAC200) 
K = CN=-NRINT—-1)/2 


SEARCH FOR KNOT INTERVAL TCNUMBER#K) C= X C= TCNUMBER*K 41) WHERE 
FPINTCNUMBER) IS MAXIMAL ON THE CONDITION THAT NROATACNUMBER) 


NOT EQUALS ZERO. 
FPMAX = QO. 
JBEGIN = 1 
DO 20 J=1i.NRINT 
~ JPOINT = NRDATACJ) 


IFCFPMAX.GE.FPINTCJ) «OR. JPOINT.EQ.0) GO TO 19 


FPMAX = FPINTCJ) 
NUMBER = J 

MAXPT = JPOINT 
MAXBEG = J3EGIN 


TZ 


5790 
5800 
po U8 
5820 
2630 
5840 
5850 
5860 
5870 
5880 
5890 
3900 
5910 
5920 
5730 
5940 
5950 
5960 
5370 
5980 
5390 
6000 
6010 
6020 
6030 
6040 
4050 
6060 
6070 
$080 
6090 
6100 
6110 
6120 
6130 
6140 
6150 
6160 
6170 
6180 
6190 
6200 
6210 
6220 
6230 
6240 
6250 


6270 
6280 
6290 
6300 
6310 
6320 
6330 
6340 
6350 
6360 
6370 
6380 
43290 


10 JBEGIN = JBEGIN¢JPOINT 1 
20 CONTINUE 
C LET COINCICE THe NEw KNOT TCNUMBER*K*1l) WITH A DATA PSINT XCNRX) 
C INSIOE THE OLO KNOT INTERVAL TCNUMBER*K) C= X C= TCNUNBEReKRO1). 
IHALF = MAXPT/241 
NRX = MAXBEG¢IHALF 
NEXT = NUMBER] 
IFCNEXT.ZGT.NRINT) GO TO 40 
C ADJUSTS THE OIFFERENT PARAMETERS. 
D0 30 J=NEXT,NRINT 
JJ = NEXT*NRINT=J 
FPINTCJJ*1) = FPINTC JJ) 
NROATACJJ¢1) = NROATAC IJ) 
JK = JSR 
TCO JR*1) = TCIK) 
30 CONTINUE 
40 NROATACNUMBER) = I4ALF-1 
NRDATACNEXT) = MAXPT-IHALF 
FPINTCNUMBCR) = FPMAXSFLOATCNROATACNUMBER) JSFLOATCMAXPT ) 
FPINTCNEXT) = FPMAXSFLOATCNRDATACNEXT) )ZFLOATCMAXPT ) 
JK = NEXT*K 
TCJK) = XONRX) 
N = Ne] 
NRINT = NRINT#1 
RETURN 
ENO 
SUSROUTINE OISCOCTN»K2,B) 


C SUBROUTINE OISCO CALCULATES THE OISCONTINUITY JUMPS OF THE KTH 
MER LVATIVE OF THE B-SPLINES OF OEGREE K AT THE KNOTS TCK*2)..TOCN-K—-1) 
C THE FIRST OIMENSION SPECIFICATION OF THE MATRIX EB MUST GE THE SAME AS 
G IN TRE CALLING PROGRAM? H MUST HAVE A OIMENSION SPECIFICATION AT 
Mee EAST 24k e2. 

OIMENSION TCN) 486290 .K2).nHC12) 

Kl = K2-1 

K = Kl-1 


NK1 = N=-K1 

00 60 L=K2,NK1 
EMK = €=f1 
90 10 J=1,K1 


IK = J¢Kl 
LJ = LJ 
LK = LJI-K2 


MOI) = TCLI=TCLK) 
HCIK) = TCLI-TCLJ) 
10 CONTINUE 
LP = LMK 
00 30 JjJ=1eK2 
JK = Jer 
PROG = 1’. 
00 20 I=JeJK 
PROO = PROOSHC(I) 
20 CONTINUE 
LK = LP¢K1 
BCLMK,J) = CTCLKJ-TCLP))/PROO 


LP = LPel 
30 CONTINUE 
£0 CONTINUE 
RETURN 


ENO 
FUNCTION RATIONCP1,F1,P25F25P3,F3) 
Reo venN THREE POINTS CP1,Fl).¢€F2,F2) AND €P3.F3),. FUNCTION RATION 


uaz) 


6490 
6410 
6420 
6430 
6440 
6450 
64450 
6470 
5430 
6490 
6500 
6510 
6520 
6530 
$540 
6550 
6560 
6570 
6530 
0590 
6600 
6610 
6620 
66 30 
6640 
6650 
5660 
6670 
6630 
6690 
6700 
6710 
4720 
6730 
6740 
6750 
6760 
6770 
6730 
6790 
6890 
6810 
6820 
6830 
6840 
6350 
6860 
6870 
6880 
6890 
6900 
6910 
6920 
6930 
6940 
6950 
6950 
6970 
6989 
6990 
7070 


an 


AQAAANAAAAAAAAAAANAANANAANANAAAANAAAN 


GIVES THE VALUE OF P SUCH THAT THE RATIONAC INTERPORATINGer UNG een 
OF THE FORM RCP) = CURPeVI/SCP4W) EQUALS ZERO AT Pe 
IFCP3.GT.0.) GOD TO 10 
VALUE OF P IN CASE P3 = INFINITY. 
P= (PISCFI-F3)SF2—-P2aCR2—-F 3 Der TO 7CCEIL—E 2 SF 30 
GO TO 20 
VALUE OF P IN CASE P3 “= INFINITY. 


LO Hie Sate la C c= 3) 
HZ = F2ZSCF3-Fl) 
H3 = F3%CF1-F2) 
P = ~-(€P12P2¢H34P24P3eH1 +P 32P1SH2)/CPLEH1 +P2%H2 +P 34H 3) 


AQOJUST THE VALUE OF PisF1lsP3 ANDO F3 SUCH THAT Fi >) 0 ano 3a 


20 IFCF2-LT.0-) GO TO 30 


PS pee 
Fil = F2 
GO TO 40 


30 \ Parad 


F3 = F2 


40 RATION = P 


RETURN 

ENO 

FUNCTION OERIVCTsNeCeoNKlsNUgARGoL) 
FUNCTION OCERIV EVALUATES A SPLINE SCX) OF OEGREE K WHICH IS 
GIVEN IN ITS NORMALIZED B-SPLINE REPRESENTATION OR CALCULATES 
QOERIVATIVES OF ANY SPECIFIED OROER NU. 


CALLING SEQUENCE 
VALUE = DERIVCTSN,CeNK1SNUSARGoL) 


INPUT PARAMETERS: 
T > ARRAY,sMININUM LENGTH Ny» WHICH CONTAINS THE POSITION 
OF THE KNOTS OF SCX)eI-E~ THE POSITION OF THE INTERIOR 
KNOTS TOCK¢2) eee. TON=K-1) AS WELL AS THE POSITION OF THE 
KNOTS TCl) sees TCK*+1) ANDO TOCN-K),-2-TCN) WHICH ARE NEECEDO 
FOR THE B=SPLINE REPRESENTATION. 


N s INTEGER VALUc GIVING THE TOTAL NUMBER OF KNOTS OF SCX). 

Cc > ARRAY,LENGTH NK1, CONTAINING THE B-SOLINE COEFFICIENTS. 

NK1 s INTEGER VALUE,GIVING THE OQIMENSION OF SCX),I-E-NK1 = N-K-le 
NU * INTEGER VALUE WHICH GIVES THe ORCER OF THE OERIVATIVE. 

ARG +: REAL VALUESsGIVING THE VALUE OF THE ARGUMENT. 

L >: INTEGER VALUE WHICH SPECIFIES THe POSITION OF THE ARGUMENT 


lee. TCL) <= ARG € TCLe1) OR 
L = NA1 IF ARG = TCNK141). 


OUTPUT PARAMETER: 


VALUES REAL VALUE.SGIVING THE VALUE OF THE NUTH DERIVATIVE OF 
SCX) AT K = ARG. 


OTHER SUBROUTINES REQUIREO: NONE. 


RESTRICTIONS: 
1) NU >= 0 
2) TCK*1) <= ARG C= TCNK1¢1) 
THE OIMENSION SPECIFICATION OF THE ARRAY H MUST BE AT LEAST Kel. 
OIMENSION TOCN),COCNK1) HCE) 
OERIV = 9. 
Kl = N-NK1] 
IFCNU.LT.O .OR. NULGE-K1) RETURN 
DO 100 T=1,X1 
IK = L+{i—Ki 
HCI) = CCERD 


AyeAe: 


7010 
7020 
7030 
7040 
7050 
7060 
T1070 
7080 
7090 
7100 
7110 
7120 
T1306 
7140 
TSO 
7140 
7170 
7180 
7190 
7200 
7210 
7220 
7230 
7240 
7250 
7250 
7270 
7280 
7290 
7300 
T3210 
7320 
7330 
7340 
12560 
1350 
7370 
7380 
7390 
7400 
7410 
7420 
7430 
7440 
7450 
7460 
7470 
7480 
7490 
T3500 
Tsao 
13920 
7530 
7540 
1550 
7560 
T3aTU 
7580 
7590 
T7609 


x%eaAN 


100 


200 


300 


400 
500 


600 


CONTINUE 
BREGNUSEO. 0) GO TO. 300 
NUL = NUel 
OC 200 J=2,NU1 
00 200 JJ=J,K1 
I = JSJ*#R1 JJ 
LI = L*I-Kl 
LJ = LelI-Jel 


Melee ecrC I) =HCI=19D7 CTC es)-1T CLI 


CONTINUE 
BPONUsEG.Ki-1) GO 10 500 
NU2 = NU+¢2 
OC 400 J=NU2,K1 
DO 400 JJ=JeKl 
I = J*Kl-JJ 
LI = LeI-Ki 
CJ = Lel-J+i 


ene CARG-> TCL II ISeHCI Ie CTCLID-ARG)sHCI-1))/CTCLID=TCLI)D >) 


CONTINUE 
CERIV = HC(A1) 
IFCNULEC.0) RETURN 
00 600 I[=1,NU 
DERIV = DERIV*FLOATCK1-I) 
CONTINUE 
RETURN 
ENO 


ae, 


7620 
7630 
7640 
7650 
T7660 
7670 
7680 
7690 
7700 
10 
Alera 
7730 
7740 
7750 
7760 
7770 
7780 
7730 
7800 
7810 
7820 
78390 
7840 
7850 
7860 
7870 


APPEND ES») 
THE PROGRAM TO FIND B-SPLINE COnPEL Gras 


10-Oec-13986 16:37:51 
Source Listing 10-O0ec-1984% 16:37:39 


PROGRAM 3SPLINECINPUT,OUTPUTSINFILE,OUTFILE)S 


TYPE 
me = 0..2555 
IMAGEROW1 = PACKEO ARRAY C0..255) OF BYTE; 
ROW1 = PACKEO ARRAY €0..257) OF BYTE; 
Beare = PACKEO ARRAY £€C0..-53,0..255]) OF BYTE: 
OAL = PACKEO ARRAY C0..5123 OF REAL? 
OA2 = PACKEO ARRAY (1.23003 OF REALS 

VAR 
ReCeFOeFleF2,F3eF4,F5°6,FT,FB Ps BYTe: 
F 2: ARRAY €0..65]) OF ROW]: 
AeIMAG= : SHIP; 
INFILE =: FILE OF IMAGEROW]: 
SPXeSPY,SPX1,SPY1 3 PACKED ARRAY C0..512) OF REALS 
RIeCL»AR2,C2eCOIRelLeSe JL eCOUNT,NEG SINTEGER:; 
TEMPX,TEMPY,RA = REAL; 
Me ITOPTsK sIPARgNeg ITER ep ANS» NU, ANSI GANSZ2,ANS3 *& INTEGER; 
NRLZNENOL MET SINTEGERS 
WeZ :sPACKEO ARRAY (€0..2£512) OF REAL; 
S$eZB,ZEeFP = REALS 
ARG» THETA,TOLE = REALS 
FLAG_BEGIN,FLAG_ENOeAK,LUMP © INTEGERS 
BEGsEN 3 PACKED ARRAY €1.25) OF REAL: 
BEGN,ENN = PACKEO &RRAY £1.25) OF INTEGERS 
TeCXeCY sPACKEO ARRAY €1.2-300) OF REAL: 
COL»eROW :PACKEO ARRAY £0..512) OF REAL; 
Prot) sPACKEO ARRAY C0.-512) OF REAL: 
OCY»eOCX,ARE,CY1,O0MAX,0OMIN = REALS 
CYMIN,CYMAX,MAXKCY =: PACKEO ARRAY (1.2.5) OF REAL: 
AREA =: PACKEO ARRAY €1.2£5] OF REAL? 
OUTFILE =: FILE OF IMAGEROWL: 
NAME 3 PACKEO ARRAY £€1.-20]) OF CHAR; 
PEAK,STAR,TER :PACKED ARRAY (€1..5) OF REAL; 


Cz FILTER THE POINTS #&) 


PROCSOURE STORE: 
BEGIN 
TEMPXS=C-1s TEMPY2=64-(R-1)5 
IF I>0 THEN SEGIN 
FOR jJ?=0 TO M OO BEGIN 
IFCCOLCJSI=TEMPX) ANDO CROWCIJI=TEMPY) 


THEN L2=J; 

ENOS 
COLCIJ s=TEMPX; 
ROWCIJ s=TEMPY; 
Ms=I¢ 
IT:slel; 

ENO 

ELSE BEGIN 


COLCIJ s=TEMPX: 
ROWCI}) s=TEMPY; 
M2=0;5 
Is=l; 
ENO; 
ENO; 


es 
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PROC 
BEGI 
FO 


F3IS=FCR41,C #195 


F6 
CA 


yource Lasting 


EQURE CMOVES 

N 
s=FCR-1,C]s FIs=FCR-1,C*1 ):s F 
FLS=FCR41,CI93% 

s=FECRCH1])3 FIS=FCR-1,C-1])3 

SE COIR OF 

0: BEGIN 


10-O0ec-138¢ 16:237:51 
10-Dec-1984 16:37:39 


s=FCR,C#l)$ 
F5S=FCR41,¢C—-1];3 


IF €F3=0) AND CF2=255) THEN BEGIN C2=C415 


COIR: =135 
ENO 
ELSE IF CF2=0) AND CF1=255) THEN 
STORE: C2=C*#13 COIR2 203: END 
EL.SE IF CF1=0) AND CFO0=255) THEN 
COIR?3=0- END 
ELSE IF CFO=0) AND CFT=255) THEN 
STORES R2=R-1s CDOIR22=33 ENG 
ELSE IF CFT=0) ANG CF6=255) THEN 
COIR:=33 ENO 
ELSE IF «€F620) ANG CF54255) THEN 
STORE. s=C-13; COIR>2=33 ENDO 
ELSE BEGIN ReeRels COIRS=23 ENDS 
ENO; 
Is BEGIN 


BEGIN 


BEGIN 


BEGIN 


BEGIN 


BEGIN 


IF CF5=0) ANDO CFS=255) THEN BEGIN R2=R¥1;5 


COIRs=23. END 
ELSE IF CF4=0) AND CF3=255) THEN 
STORES R2=Reils COIR =13 ENO 
ELSE IF CF3=0) ANG CF2=255) THEN 
COIRs=13: END 
ELSE IF CF2=0) AND (€F1=255) THEN 
STORE; s=Ce1$s CDIR2=03 END 
ELSE IF CF1=0) AND CFO0=255) THEN 
COIR:=03 ENO 
ELSE IF CFO=0) AND CFT=H=255) THEN 
STORES Ri=R-1:s COIR2=0% END 
ELSE BEGIN Cs=C-1ls COIR?2=33 ENO; 
ENOS 
2: BEGIN 


BEGIN 
BEGIN 
BEGIN 
BEGIN 


BEGIN 


Cs=C+l; 
Cs=C+15 
RS=zR-1ls 
Ro=eR—-ls 


C:=C-13 


IF CFT=0) ANDO (F6=255) THEN BEGIN C2=C-1;5 


COIR:=33; ENO 
ELSE IF CF6&=0) ANG CF5=255) THEN 
STORES Cs=C-13 COIR: =33 END 
ELSE IF CF5=0) AND CF4=255) THEN 
COIRS=2:; ENO 
ELSE IF CF4S=2=0) ANDO CF3=255) THEN 
STORES R2=Re13 COIRS=13 END 
ELSE IF CF3=0) ANG CF2=255) THEN 
COIRs=13s ENO 
ELSE IF CF2=0) ANDO CF1=255) THEN 
STORE? Cs=Ce1s COIR2=13 END 
ELSE BEGIN R:=R-13 CDIR:2=03 ENOS 
ENDS 
3: BEGIN 


BEGIN 


BEGIN 


BEGIN 


BEGIN 


BEGIN 


Rs=R+4+15 
Rs=Rel1l;s 
Cs=€41; 
Cs=C+15 


RS=R—-1;5 


IF C€F1=0) ANDO CFO=255) THEN BEGIN R2=R-1, 


COIRS=0. ENO 


ELSE IF CFO=0) ANC CFT=255) THEN BEGIN C2=CH-l1. 


TZ 


VAt=T 
_~ORAOD 
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STGREs R2=R-1s COIR:S=0%3 ENO 

ELSE IF CF7=0) AND CF6=255) THEN BEGIN C:=C-1; 
COIR:=33; ENO 

ELSE IF CF6=0) ANDO CFS=255) THEN BEGIN R2=R415 
STORES s=C-1; COIR:=2: ENO 

ELSE IF CFS=0) ANO CF4S2=255) THEN BEGIN R2=R41;5 
COIR:=2: ENO 

ELSE IF CF4=0) ANO CF3=255) THEN BEGIN C2=C+15 
STORE; s=R+1; COLIR2=2. ENO 

ELSE BEGIN Cs=Cel1s COIR: =1;. ENO: 

ENO; 
ENO; 
END; 


PROCSQURE INITIALS 


BEGIN 
FOR I :=0 TO 255 OO BEGIN 
XCIJQ :=0; 
YCIZ :=0; 
ENO; 
I:=0; 
ENO; 


PROCEDURE FIRSTS 


VAR 
NeoM sINTEGERS 
SEGIN 
N 3= Of M2=05 
FOR C :=0 TO 200 OO SEGIN 
FOR R 3:=0 TO 63 OO BEGIN 
IF CFOR,CI=255) ANDO CN=0) THEN BEGIN 
Ci s= Cs RI 2=R3 N S= Nels 
END: 
IF CFOR,.255-CIJ=255) AND (M=0) THEN BEGIN 
C2 3:2 255-Cs$ R2 2= RS M S=MOl§ 
ENO S$ 
ENO; 
ENO; 
ENO: 
PROCEDURE ROTATES 
BEGIN 
NEG:=0; 


IF R1=R2 THEN BEGIN THETA:=0.03 
FOR I:=0 TO m OO BEGIN 
XCIJ°=COLCIIJ;s YCOIJ2:=ROWCII:S 
ENO; 
ENO 
ELSE THETASSABSCARCTANCCR2—R1)7¢€CC2-Cl1)))3 
IF CTHETACIO) ANDO CR2>R1) THEN BEGIN 
FOR I:=0 TO ™ OO BEGIN 
XCIJS=COLCIJSCOSCTHETA)—-ROWCIISSINCTHETA): 
TET IS=ECOLELISSINCTHETAD+ROWLIDSCOSCTHETA): 
IF YCI3>=63.0 THEN NEG3=NEG412 
ENO; 
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END 
ELSE IF CTHETACS0O) ANDO CR2CR1) THEN SEGIN 
FOR I:=0 TO ™ OO BEGIN 
XCIJS=COLCIISCOSCTHETA)*ROWCIIJSSINCTHETA), 
YCOIJS=-COLCIISSINCTHETA) +ROWLCIJ=COSCTHETA): 
IF YOII>=63.0 THEN NEGS=NEGC1; 
ENO; 
ENOS 
C* FILTER *#) 
I:=0; 
FOR K:=0 TO M OO BEGIN 
TEMPXS=XCKIs TEMPYS=YCKI; 
IF I>0 THEN BEGIN 
FOR J:=0 TO M OO BEGIN 


IF CXCJJ=TEMPX) AND CYCIJ=TEMPY) THEN iieJ3 

END: 
XCIJQS=TEMPXs YCIJS=TEMPYs Meals ITs=elel;s 

ENO 

ELSE BEGIN XCIIS=TEMPXs YCIJS=TEMPYs MF=03 Ti=13 

ENOS 

ENO; 
ENO; 


PROCEQURE PARAMC X20A1s Y2OA1l;s VAR Z20A1s Wi0Al; 
MSINTEGERs VAR ZBZREALS VAR ZESREALS KFINTEGERS 
SsREALS VAR NSINTEGERS WAR T20A23s VAR CX:20A2;5 
VAR CY2OA23 VAR FPSREALS IOPTSINTEGERS 
IPARSINTEGER?S VAR TERSINTEGER)$ FORTRAN, 

PROCEQURE INITTCSPEEOQZS INTEGER): FORTRAN; 

PROCEQURE VWWINOOCXMINS REALS XRANGESREALsS YMIN: REAL; 
YRANGESREAL)S FORTRANS 

PROCEQURE MOVEACXSREALS YSREAL)D: FORTRANS 

PROCEOQURE ORAWACXSREALS YIREAL)DS FORTRAN, 

PROCEQURE ANCHOCICHARZ INTEGER): FORTRAN; 

PROCSOURE FINITTCILS INTEGERS I22:INTEGER): FORTRAN, 

PROCEQURE OASHACKSREALS Y2REALs LSINTEGER): FORTRAN;G 


FUNCTION OERIVCZREF 20822 NZINIEGER. @ CX DA, 
NKLSINTEGERs NUZINTEGERS ARG:SREAL; 
LSINTEGER) REALS FORTRAN: 
PROCEQURE SPLCOEFs: 
BEGIN 
7=53 LUMP:=0; 
WHILE CI<=N-5) ANO CLUMP=0) OO BEGIN 
IFCCYCI-1LIQ>CYCTII) ANDO CCYCT*¢13>9Cr¥C1I1) AND 
CCYCTe¢22>9CYC1¢13) THEN 
LUMP 3=13 
Ts=I¢ls; 
ENO; 
IF LUMP=1 THEN BEGIN 
FLAG_BEGIN:=03 AK2=13 FLAG_ENOS=03 [:=53 
WHILE I<¢=N-5 OO BEGIN 
IF FLAG_BEGIN=0 THEN BEGIN 
IF CCYCI-13>CYCII) ANO CCYCI*1I359CYCII) ANDO 
CCYCI*¢ZI5CYCI 13) THEN BEGIN 
SEGCAKJS=TCII: FLAG _SEGINi=135 SEGNCAKJ2=1; 
ENO. 
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END 
ELSE IF FLAG BEGIN=1 THEN BEGIN 
TF C€VEI=$19>=CYEI3) ANDO CCYECIIJ>=CVET*13) AND 
CCYCI¢23>=CYCI4¢13) THEN BEGIN 
ENCAKJ2=TCI+4139$ FLAG_ENOs=1s ENNCAKJS=Tels 
ENDS 
ENDS; 
IF CFLAG BEGIN=1) AND CFLAG_END=1) THEN GEGIN 
FLAG _ BEGIN: =03s FLAG_LENDS=05 AKS=AKe1s Li=I-1; 
ENDS 
Is=Iels; 
ENDS 
END 
ELSE IF LUMP=0 THEN BEGIN 
FLAG_BEGIN:=0; FLAG_END2=0- I2=S5s AK2=1;5 
WHILE I<=N-5 OO BEGIN 
IF FLAG BEGIN=0 THEN BcGIN 
TECEYE T-175>=CYEiI3) AND CCYEI 41 2>CYCID) And 
CCYEL*2Z)>CYE1I]) THEN BEGIN 
BEGCAKJ2=TCIIJ: FLAG_BEGIN: =13 BEGNCAKJ:=I; 
ENDS 
SND 
ELSE IF FLAG _ BEGIN=1 THEN SEGIN 
TEC EYE I=129=CYCII) AND CCCYCID<=CYCI¢13+TOLe) 
AND CCYCIJD=CYCI*¢1I3I-TOLE)) AND 
CCCYCII<=CYCI+¢23+TOLE) AND 
CCYCIIJD=ECYCI*+23+TOLE)) THEN’ BEGIN 
ENCAKJS=TCIIJ:s FLAG_ENDS=1:3 ENNCAKJS=I5 
END: 
END; 
IF CFLAG_BEGIN=1) AND CFLAG_ENOD=A1) THEN JEGIN 
FLAG_BEGIN:=03 FLAG_END2=03 AKS=AKRe13 [221-13 
ENG; 
PeH=t ers 
ENDS 
IF FLAG_BEGIN=1 THEN BEGIN 
I: =BEGNC AK]; 
WHILE CIC¢=N—-5) AND CFLAG_END=0) OD BEGIN 
Perce YEl—1J>=CYCII) AND 
CCCYCII<=CYCI+¢1I3+TOLE) 
AND CCYCIJD=CYCI*1I-TOLEI) THEN BEGIN 
ENCAKJS=TOCIINs ENNCAKIJS=Is FLAG_END2=15 
ENO; 
Ts=1413% 
ENDS: 
END; 
END; 
ENDS 


PROCEDURE ALUMP? 
BEGIN 
Ake =) * 
WHILE CENNCAKJ)>O DO BEGIN 
AREATAKIS=03 CYMAXCAKI2=0.93 CYMINCAKI2=1000.03 
FOR JT2=CBEGNCAKRI) TO CENNCAKJ) OD BEGIN 
IF CYMAXCAK) <= CYCIIJ THEN BEGIN : 
CYMAXCAKJS=CYCIIDs MAXCYCAKJS=STCII 
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END; 
IF CYMINCAKD d= CYLII THEN 
CYMINCAKJ:3=CYCII3 
ENO; 
CY1L2=CYCBEGNCAKJI-CYMINCAK IS 
FOR L2=CBEGNCAKIQ) TO CENNCAKI-1) OO BEGIN 
OCYsscCYCI+#13-CyYcC1lI; 
VEXS=CXEL lel IJ-CxLl as. 
AREALCAKJ2=AREACAKI+CY1L20CX*OCYsOCx/2.03 
WRITELNC°AREA=",AREACAKI,” J="oIy” AK=",AK)3 
CY¥13=CvY1+0CY3 
CNO3 
AK S=AK413 
ENOS 
ENOS 
PROCEOQURE POSINX$ 
BEGIN 

13203 MET:=03 

WHILE CI<=M-1) AND CMET=0) OO BEGIN 
IF TeEmMPx=ZCI3 THEN BEGIN 

TEMPYcS=XCIIs METS=A15 
END; 
s=Iels 
ENDO; 
ENOS 
C£203452 342255252 -- 45a oe eo ee 
Ct MAIN PROGRAM 2) 
BEGIN 

WRITELNC “INPUT CUT OR CONLINE FILENAME” )3 

REAOLNCNAME); 

OPEN CINFILEsNAME,HISTORY = OLD, 
ACCESS_METHOO :s=SEQUENTIAL, 
RECORDC_LENGTH :=256,RECORO_TYPE s=FIXED); 

RESETCINFILE)? 

R segs 

WHILE NOT EOF CINFILE) OO 

BEGIN 
READ CINFILEsIMAGECRI): 

FOR C := 0 TO 255 OO 
FOR+#1,C+#1) s= IMAGECR,CI: 
RQ s= Rei 

ENO; 

CEOS eCiN@ Tee). 

FIRST: 

R S#R1¢ °C 2=C€1i; COIR 221; 

INITLA&L: 

WHILE C€C<¢>C2) OO BEGIN 
STORE: 

CMOVES 

ENO$ 

ROTATE: 

FOR 13:=0 TO ™ OO BEGIN 
IF NEGCSO THEN YC 72=7679-20_0. 

END; 

MS=Mo1s KSe2e SP=O0518 > TORT? =0% ean 0. 

FOR I:=0 TQ m™ OC 
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Aes = 150% 
ANSs=13 
WHILE ANS=1 O00 BEGIN 
WRITELNC°OLO S=",S,°NEW S=°)3 
READCS): 
WRITELNC“°OLO K=79Ky"K=")3 
REAQCK): 
PARAMC Xe VapZaoWeoMepZBeZeEnKe So NygTeCXeoCYeFP,IOPT, 
IPAR.IER); 
MARITELNC® S=°,S5e” IER=°, ITER,” Ma" eM,” NW=°QN, 
° CXCNI=",CXONI, ° CYON-439=°,CYON-4))3 
WRITELNC” PRINTING YVES=1 PLOT X-Y=2,CX,CY=3 °, 
“X~Y-CX-CV¥=6°,° NO=5°)3 
REAQCANS1)3 
IF ANS1=1 THEN S8EGIN 
FOR I2=1 10 N OO 
MRT EIUNG Che oe CXL a. CY¥="“.CYET; 
= T=*, 701293 
ENO 
ELSE IF ANS1=2 THEN BEGIN 
PRITIC9I60O); 
VWINOOCO.0,20¥"-13,0.0,256.0)3 
MOVEACZC03,XC0]); 
FOR I:=0 TO M=1 QO 
ORAWACZCIIJACIID: 
MOVEACZCO0I3-YCO0I): 
FOR I:=0 TO M=-1 OO 
GASHACZE PI YoLI,2)$ 
FINITTCO,0)3 
ENO 
ELSE IF ANS1=6 THEN SEGIN 
INITTC960)3 
VWINOOCO0.0.Z0M-13,0.0,256.0)3 
MOVEACZCOI,X003);3 
FOR 12:20 TO m=1 OO 
ORAWACZCII,X£1I);3 
MOVEACZCOI,YL0I)3 
FOR I:=0 TO M-1 OO 
ORAWACZCII,YCII): 
MOVEACTOC4),CKO49)3 
FOR IT:=4 TO N-s& OO 
GOASHACTOCII.CXLII3.,2)3 
FOR *=4 TO N=4 OO BEGIN 
POVERQTETI=-1.5;,CXCII3-1.5)5 
ANCHOC111)3 
ENCS 
MOVEACTESI,CYC4I3)5 
FOR I:=4 TO N-4 OO 
DashHactTCla,cYCI +295 
FOR I:=46 TO N=4 OO BEGIN 
MOVEACTCII-1.59,CYCII-1.5)3 
ANCHOC111):3 
ENO; 
OMIN: 20.03 OMAX2=256.0: TEMPX:2=OMIN? 
WHILE TEMPX<=ZCM-1)3 OC SEGIN 
MOVEACTEMPX,OMIN)' 
ORAWACTEMFX,OMIN)D; 


oe 


10-Dec-1986 16:37:51 VAX=] 
Source Listing 10-Vec-198% 16:37:39 _ORAC 


DRAWACTEMPX,DMAX)> 

TEMPXS=TEMPX4+10.03 

END; 

OMIN2=0.03 OMAX2=ZCM—-1)5 TEMPX>S=DMING 

WHILE TEMPX¢€=256.0 OO SEGIN 
MDVEACDMIN,TEMPX); 
DRAWACOMIN,TEMPX)>: 
ORAWACDMAX,TEMPX)S 
TEMPXS=TEMPX410.03 

END; 

FINIT TCO. Co. 

ENO 
ELSE IF ANS1=3 THEN BEGIN 

INITTC960);5 

VWINDDCO.0,TENI,0.02256.0)3 

MOVEACTEC4),Cx04))5 

FOR [:=4 TO N-4& DD 
ORAWACTCIQ,CXCIID): 

FOR I:=4 TO N~-S OO BEGIN 
MOVEACTEILIJ=-1 2 5,CXL1I-1.5) 
ANCHOC111)3 

END; 

MOVEACTC4SI,CYC4)); 

FDR I:=4 TD N~-4& DO 
OASHACTCII.CYCII.,2); 

FDR I:=4 TD N—-46& OO BEGIN 
MOVEACTECII-1.5,CYCiIjJ-1.5)35 
ANCHOC111)3 

ENO; 

FINITTCO,0)3 

ENO; 

Cx IMPORTANT FDRTRAN DECLEAR FROM 1 *) 
NKITS=N-K-13 Lo=Kels NENOSH=N—-KS J12=05 
FOR s=L TO NENO OO BEGIN 

ARG: =TCTII: 

WHILECARGO=TC L419) AND CLCNK1) DO 

s=tLe1, 

SPK1ICJ1I:2:= OERIVCTN»eCX»NK1,0.ARG OLDS 

SPY1ITILIS= OERIVCTNeCY eNK1,0,4RGeL)S 

Js =) 1425, 

ENO; 

J2=O05 Li=kKels 

FOR s=L TD NENO CO BEGIN 
ARG:=TCI2;3 

WHILE CARG>=TEICL +413) ANDO CLCNK1) OO 

Ls=Lel; 

TEMPXS=TCI*¢1I~AaRGs 

IF TEMPX>=7.0 THEN BEGIN 
WHILE ARG<TCI4¢13 OO BEGIN 

SPXCJJ 2= DERI VC Tee CXeNK1,0,ARGoL)S 
SPYCJI = DERIVCTeNeCY,NKR1eO0eARGot)s 
ARG =ARG4*2.05 Js=Jels 


END. 

END 

ELSE SEGIN 
SPXCJJ 2= DERT VOT eo CX eNK1,0,ARG,L); 
SPYCJIJ = DERIVCTNeCY»NK1,0,ARGoL); 
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8 Severs |e Ts 
ENO; 
ENG; 
WRITELNC°O0 YOU WANT PLOTTING YES=1 “es 
"“PRINT=2 COMPARE=3 NOQ=4");5 
REAQCANS2)>; 
IF ANS2=1 THEN BEGIN 

INITTC960)3 

VWINODO0C0.02256.0.0.0+256.0)3 

MOVEaCXEGS.,YLO3); 

FOR I1:=0 TO M=1 OO 
ORAWACXCIIJ,YCIID: 

FOR [:=0 TO Jl-1 OO BEGIN 
MOVEACSPRICID—-125.SP41CII-1.5); 
ANCHOC111)¢ 

cNO$ 

MOVEACSPXC0],SPYCTOI); 

FOR I:=1 TO J=-1 OC 
CASHACSPXCIIJ.-SPYCII.2):5 

SINITITTCOs0)> 

END 
ELSE IF ANS2=2 THEN SEGIN 
FOR I:=C TO Jl=-i OO 
WRITELNC°SPX1=°,SPX1CII, 
. SPY1=" »SPY1CII)3 
=NO 
ELSE IF ANS2=3 THEN BEGIN 

WROEEELNC"TOL=" 25 

REAODCTOLE): 

SPLCGEF¢ 

ALUMP 

AKS=13 RAS=XCM—-1I-xXC0I3 

WHILE ENC AKIJ>DO OO BEGIN 
TEMPXS=S=BEGCAKIJ$3 POSINX: 
STARCAKJS=CCTEMPY=XC0) J—-RA/S2)/RAS 
TEMPX2=ENC AKI: POSINX: 
TERCAKJS=CCTEMPY=XCOJI-RASZ)/SRAS 
TEMPX2=MAXCYCAKI: POSINX: 
PEAKCAKJ2=CCTEMPY—-XC03)J—-RA/S2)/RAE 
AREACAK] sS=CAREALAKI)/CRASE2 D3 
WRITELNC°BEGIN=E"STARCAKI,” ENO=", 

TERTCAKI,” TOTAL=",RA, 
" AREA=",AREACAK), 
- PEAK=",PEAKCTAKI))3 
AKS SAK +15 
ENOs 
ENO: 
WRITELNC "O00 YOU WANT RUN AGAIN YES=1", 
S enNG =< js 
REACCANS): 
IF ANS=1 THEN EBEGIN 

WRITELNC"°IOPT=")35 

REACCIOPT): 

ENO; 
ENO$s (€* WHILE =) 
ENG. 
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